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Abstract. This paper develops results for the next nearest neighbour Ising model 
on random graphs. Besides being an essential ingredient in classic models for 
frustrated systems, second neighbour interactions interactions arise naturally in several 
applications such as the colour diversity problem and graphical games. We demonstrate 
ensembles of random graphs, including regular connectivity graphs, that have a 
periodic variation of free energy, with either the ratio of nearest to next nearest 
couplings, or the mean number of nearest neighbours. When the coupling ratio is 
integer paramagnetic phases can be found at zero temperature. This is shown to be 
related to the locked or unlocked nature of the interactions. For anti-ferromagnetic 
couplings, spin glass phases are demonstrated at low temperature. The interaction 
structure is formulated as a factor graph, the solution on a tree is developed. 
The replica symmetric and energetic one-step replica symmetry breaking solution is 
developed using the cavity method. We calculate within these frameworks the phase 
diagram and demonstrate the existence of dynamical transitions at zero temperature 
for cases of anti-ferromagnetic coupling on regular and inhomogeneous random graphs. 
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1. Introduction 

Many systems exhibit complicated low temperature phases that are described at a 
phenomenological level by classical spin state lattice models. A necessary ingredient 
in the description is often frustration in interactions, and one of the simplest origins 
of frustration is a competition between nearest and next nearest neighbour couplings. 
The simplest models can already demonstrate a range of phenomena with interesting 
long range correlation and modulation patterns. At the same time many problems in 
optimization have a similar structure, an objective function with solutions equivalent 
to the ground states of some next nearest neighbour Ising model. In these optimization 
problems the ground state space structure can be related to issues of algorithmic 
complexity in resolving decision problems and sampling. 

Theoretical developments based on the cavity and replica method have allowed 
greater insight into these problems in both traditional solid state physics, theoretical 
computer science and other cross disciplinary fields [U E]. This includes dilution 
effects that mimic the finite range of interaction of spins, with much recent focus on 
random graphical models. Random graph models allow dilution effects to be studied 
within a mean-field analysis. More recently, much development of the theory of the 
ground state space of Ising models on random graphs has occurred in the guise of 
constraint satisfaction problems [31 H] . Random graph models, including inhomogeneous 
connectivity, provide benchmarks for optimization problems such as kSAT and colouring, 
and a basis for powerful coding methods. 

Two famous models in physics involving frustrated couplings are the Axial Next 
Nearest Neighbour Ising (ANNNI) Model [5], and the Next Nearest Neighbour Ising 
model (NNNI) [§J. The difference between the two is that couplings in the former are 
frustrated only along a particular axis, whereas in the latter interactions are symmetric 
along all axes. The thermodynamics in the bulk depend on the ratio of the nearest to 
next nearest couplings, and the lattice structure. In these models the ground state is 
normally dominated by ferromagnetic and modulated (2) ordering. The competition of 
the nearest and next nearest interactions leads to a weakening of long range correlations, 
and hence a decreased Curie Temperature, or even the absence of a transition at some 
special coupling ratios. A simple argument based on the local balancing of ferromagnetic 
and modulated tendencies provides an accurate prediction for some of these ratios (61 [7] . 
At non-zero temperature long range correlations are changed, and other modulation 
patterns may be thermodynamically dominant, or relevant as meta-stable attractors for 
dynamics. At some special ratios one can raise the temperature and observe an order 
from disorder effect [8J. A Devil's staircase phenomenon can also be found varying the 
coupling ratio at fixed temperature [9]. 

In understanding the phase diagram, locally tree-like models have become popular. 
The Vannimenus model was proposed as the Bethe-Peierls approximation to the ANNNI 
model PHI ttU EG] , with analogous studies in the case of the NNNI model [IH1 HU EE] • 
These reproduce many thermodynamic features of the finite dimensional models 
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including special ratios. Other frustrated next nearest neighbour models have been 
studied on Bethe Lattices and Cayley trees coinciding quantitatively with ANNNI or 
NNNI models in some parameterizations. Amongst these are Husimi Cactii Ising spin 
models, Potts and continuous spin NNN models and third nearest neighbour models 

[ini nzi nn nni isni eh E2] . 

The purpose of this paper is to develop an approach to the NNNI model utilizing 
insights from hard optimization problems and the cavity method. Our original 
motivation was the diversity colouring problem, an algorithmically hard problem that 
finds its modern applications in distributed storage on networks that are typically 
disordered [231 EH EE]. A statistical physics approach for the case to the anti- 
ferromagnetic Potts model (Q-colour states) has been previously studied, demonstrating 
spin glass and paramagnetic phases [23]. In this paper the Ising spin case is studied for 
general coupling ratios, and a one-step replica symmetry breaking analysis is undertaken 
to demonstrate a more precise description of the phase space that is expected to 
generalize in many ways to Potts spin models. 

A closely related application of NNNI models is in the analysis of pure Nash 
Equilibria in graphical games [261 EZ], the distributed storage problem may also be 
posed in such a framework. In this context each node (player) on a graph may choose a 
strategy (such as ±1) so as to maximize his objective that depends on the action of his 
neighbours, this objective is often a simple pairwise function. Since the optimal strategy 
depends on the joint distribution of the neighbours' strategies there is a second neighbour 
interaction through the intermediate player. Certain types of Nash equilibria can be 
shown to correspond to the ground states of a NNN model, and other thermodynamic 
properties can have game theoretic interpretations. 

In section El we begin our analysis with a formulation of a graphical model that 
enables the solution on a tree to be presented. In section [3] we develop the free energy 
for a tree as a simple exact solution of the graphical model. In section H] we identify a 
simple relationship between coupling ratios and connectivity, identifying special ratios 
based on local properties in the graphical model. In section |5] we then present ground 
state statistics for some loopy graphs of small size sampled from the linear ensemble, 
which demonstrates clearly the significance of the special ratios. 

In section |6] replica symmetric (RS) and energetic one-step replica symmetry 
breaking (1RSB) forms of the cavity method are developed. Equivalence between our 
model and classes of constraint satisfaction problems (CSPs) are demonstrated. In 
section [7] the phase diagram is presented as a function of temperature, connectivity 
and coupling ratio |28j at the level of replica symmetry for cases of anti-ferromagnetic 
nearest and next nearest couplings. The phase diagram allows a paramagnetic phase at 
high temperature, which at the special coupling ratios can persist to zero temperature. 
At low temperature we find a spin glass solution. For graphs sampled with regular 
or linear connectivity we consider special coupling ratios and show that with variation 
of connectivity, there may exist both paramagnetic and strongly correlated extensive 
entropy phases at zero temperature. Paramagnetic behaviour is shown to follow a 
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predictable periodic pattern with mean connectivity. In regimes with a competition 
of nearest and next nearest couplings both dynamical and continuous transitions 
from paramagnetic to spin glass phases are shown. The replica symmetric and zero 
temperature one-step replica symmetry breaking solutions for the spin glass are found 
to be unstable; whereas the paramagnetic solution is the unique stable solution in some 
parts of the parameter space. We conclude with a discussion in section |HJ 

Appendices include additional numerical results on small graphs, an investigation 
of the support for the energetic 1RSB method, and details of the numerical methods 
used to solve the cavity equations. 

2. The Model 

We consider a spin state model with an interaction structure described by a graph 
G(V, E), consisting of a set V of N vertices and a set E of M edges. Each vertex % in V 
is associated with a spin variable Si = ±1. Two vertices, if connected by an edge, are 
called nearest neighbours, and denoted If two distinct vertices are both nearest 

neighbours of a third vertex, then they are next nearest neighbours, and denoted 
We consider simple graphs, so that there are no self-loops, and no double edges - vertices 
are never nearest neighbours of themselves and never twice nearest neighbours of each 
other. However, vertices may be second neighbours by multiple paths, and J\f)ij( denotes 
the number of such paths. Every path contributes additively to the coupling. In our 
asymptotic analysis we will neglect the effects of neighbourhoods in which vertex pairs 
are both nearest and next nearest neighbours (due to a loop of length three), and cases 
where N)ij( exceeds one (due to loops of length four). 

Given the set of first and second neighbours, we write the Hamiltonian 



which is an Ising spin model with the ratio of first and second neighbour couplings 
defined by A. The second neighbour couplings are always anti-ferromagnetic, positive 
A indicates anti-ferromagnetic nearest neighbour couplings and negative A indicates 
ferromagnetic nearest neighbour couplings. Frustration, the key feature of the model, 
may be present even where A = 0, though we concentrate in this paper on the case of 
positive A. 

When considering macroscopic properties we consider a graph sampled from an 
ensemble, an example is shown in figure [TJ for which we also illustrate a ground state 
assignment of the spins. We consider mainly an ensemble of random graphs subject to 
a linear connectivity distribution constraint, parameterized by mean connectivity C. A 
graph from this ensemble is equivalent to a random regular connectivity graph when 
the mean connectivity is an integer, or else contains vertices with connectivity \C\ and 
\C~\ such that the connectivity distribution is described by 




(1) 



P(C) = S Cil c } (\C] -C) + 6 c ,m{C ~[C\). 



(2) 
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Figure 1. In our model spins interact according to a random graph structure. We 
illustrate a ground state on a graph of N = 100 vertices, sampled from a linear 
connectivity ensemble. All vertices have connectivity two or three, mean connectivity 
is 2.56, and A = 1. Black indicate +1 spin assignments, and white —1. 



Here \_x\ ( \x~\ ) notation is used to denote the nearest integer smaller (respectively larger) 
than or equal to x. 

Some topological definitions are useful: a loop is a closed sequence of vertices 
connected by edges; a leaf is a vertex of connectivity one; a tree is a connected graph 
without loops. This terminology extends naturally to the factor graph defined later. 

Assuming the linear ensemble, the model can be parameterized by mean 
connectivity C, the coupling ratio A, and we also consider its dependence on the inverse 
temperature (3. 

2.1. The factor graph representation 

The factor graph representation is a bipartite graph representation of the model and can 
be derived following the cluster variational method at distance one, or equivalently a 
simple application of the region graph method [291 [30]. Alternatively, for those familiar 
with the junction tree method, one can realize that ignoring long loops the width of the 
graph is locally two - as such one can define pseudo-states to create a locally tree-like 
graph in the Markov sense - this is precisely the graph we formulate. The factor graph 
encodes all nearest and next nearest neighbour couplings amongst spins in terms of 
multi-body interactions amongst some generalized variables. Under this framework any 
neighbourhood in the graph G(V, E) without loops is mapped to a region in the factor 
graph without short loops between generalized variables; so allowing locally tree-like 
methods to be applied. 

We first transform the Hamiltonian |1| to a sum of local multi-body interactions. 
For each vertex in the graph, i, we define a multi-body interaction dependent on the 
spin Si and the set of nearest neighbour spins S nn ^y The multi-body interaction at 
factor i is 

H t (S l ,S nn{l) ) = U\S l + S j) -H?(Ci,\), (3) 

y jenn(i) J 

with nn(i) the set of all vertices that are the nearest neighbours of i, and Ci = \S nn f{\\ 
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is the vertex connectivity. The second term normalizes the local ground state energy to 
zero, 

H?(Ci,\)= min Hi(Si,S nn ®). (4) 

The Hamiltonian is 

H(S) = Y,H l (S l ,S nn{i) ) , (5) 

i 

which is equivalent to JT]) up to a quenched additive term, and is conveniently bounded 
below by zero. 

Generalized variables, which we can henceforth call g-spins, are introduced in one- 
to-one correspondence with the edges of G(V,E). A g-spin, located at edge is an 
ordered spin pair: created from copies of the spins Si and Sj. Factors are introduced in 
one-to-one correspondence with the vertices in the original graph. A factor % is connected 
by edges to the C, generalized variables {(Si,Sj)\j G nn(i)}. Factor i incorporates 
intersection constraints and the multi-body energetic term Hi ([3]). The intersection 
constraint at % requires that, for any Hamiltonian or temperature, the elements Si in 
each neighbouring g-spin are identical. Note that although the state space expands to 
(2 2 ) M in the new model, there are only 2 N assignments consistent with the intersection 
constraints as expected. Frequently we will be interested in whether the state has two 
spins aligned, which we call a dimer, or misaligned, which is called a non-dimer. Note 
the intersection constraint depends on the parity of the constituent spins, but the energy 
does not. 




Figure 2. (a) A locally tree-like subgraph of G(V,E) is demonstrated, we label 
the central vertex and its nearest neighbours; vertices are shown by circles and have 
connectivity 1, 2 or 3 in the example, (b) With circles representing spin variables 
the interactions can be represented either by two types of binary couplings (left), 
or by multi-body interactions (right). The multi-body interaction labeled i encodes 
Hi(Si, S nn (i)). Spin i interacts in three other multi-body interactions (not shown) 
since it is also in the sets S nn (j), S nn ^) and S nn (iy (c) g-spins variables can be 
associated with the edges of the graph G(V, E) and are represented as circles. Each 
factor encodes an energetic multi-body interaction, and an intersection constraint. The 
graph is locally tree-like, unlike the interaction structures implied by figure (b). 
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To summarize, the factor graph consists of a set of N factor vertices, a set of M 
variable vertices, and a set of edges representing their dependencies. Each variable 
is called a g-spin and has four states. The relationships with the original graph are 
visualized in figure [2j where factors are represented by squares and generalized variables 
by circles. The factor graph is found deterministically from G(V,E): 

1) Each factor vertex corresponds to an element in V, vertex labels i define factors 
uniquely; 

2) Each variable vertex corresponds to an element in E, directed edge labels define 
g-spin(s) uniquely 

3) An edge exists from generalized variable (j, k) to factor vertex % only if z = j or 
i — k. 

By construction every variable vertex in the factor graph is connected to exactly two 
factors, and each factor vertex is attached to Cj variable vertices. 

2.1.1. Cavity graphs The cavity graph Ti is a subgraph of the full factor graph formed 
by removing one factor i and its adjacent edges. In this subgraph g-spins previously 
attached to factor i have become exceptional: they are connected to only one factor. 
The factor graph Tj is that part of the cavity graph J~ i connected to g-spin (>Sj, 5^). 
The Ci cavity graphs derived from Ti (J~*-n) are disjoint if G(V, E) is a tree. 

To describe relationships on the factor tree we use the terminology of trees. In a 
factor graph J~j-+i we say that generalized variable vertex is the root of the graph. 
The root (ancestor) has Cj — 1 descendants, which are the vertices attached to factor 
j excluding the root itself. Factor graphs of the full model are shown in figure [2j with 
some additional notation relating to the cavity method shortly to be discussed. 

2.2. Auxiliary Hamiltonian 

In developing our analysis the variety of energy levels applicable in the multi-body 
interaction can be problematic. An auxiliary Hamiltonian is proposed and utilized in 
the 1RSB method. The local ground state energy continues to be zero energy, but all 
excitations are fixed to energy one; flH]) is mapped to the new form by 

Hi(Si, S nn (i)) = I [Hi(Si, S nn (i)) > 0) . (6) 

using the indicator function I that evaluates to 1 if the argument is true and zero 
otherwise. Note that (j3J) and (EJ) are interchangeable only if local excitations are absent. 
Thermodynamically this requires both that the temperature is zero, and that the ground 
state is of zero energy (unfrustrated) , which is not true for many of the ensembles studied 
herein. This Hamiltonian is exploited only in the zero temperature analysis within this 
paper, use of (J3J) is to be assumed unless otherwise stated. 
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Figure 3. Beginning with a subgraph of the full factor graph figure this figure 
demonstrates the construction of the free energy from cavity graph fields, (a) Knowing 
the fields on two factor graphs convergent on a single edge (i, j) these can be converged 
to give the edge free energy (TT6l) . (b) Calculation of (|T9| is possible given the Cj factor 
graphs convergent on a factor i. The probability of the root state marginalised over 
all other g-spins is encoded in the vector in a tree. 



2.3. Cut-Poissonian connectivity distribution 

The Cut-Poisson ensemble is also considered in experiments. The marginal probability 
distribution is given by 

P(C)^(C>2) ^ 2 - g _>g- 2 » C , (7) 

The connectivity at any site is two plus a Poisson distributed random number. 

Both ensembles presented in (jTJ) and ([2]) have minimum node connectivity two, and 
have modest variance in node connectivity. The ensembles show a degree of variability 
in results, whilst avoiding some more simple sources of variation such as the size of graph 
cores. Henceforth use of the linear ensemble (j2J) is to be assumed unless otherwise stated. 



3. Free energy on a tree 

The partition function where G(V,E) is a tree can be determined recursively and 
efficiently by the Bethe-Peierls recursion. The cavity method allows an extension of 
this method to locally tree-like graphs in section |6j This formalism is equivalent to that 
employed in [23], in which it was applied to Potts spins with A = 1. 

The partition function can be factorized about an arbitrary g-spin (S^, Sj) , using 
g-spin probabilities and partition functions on cavity trees. The two cavity graphs J~i-+j 
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and J-j-+i are trees, and non-interacting given (Si,Sj). Thus the partition function of 
the full graph is 

Z(T) = Z^j^T) = 22 [Zi->jPi->j{Si, Sj)) [Zj->iPj-+i(Sj, Si)) , (8) 



Z^j is the partition function on the corresponding cavity tree, summing over all g-spins 
including (S i7 Sj). We define Pi^j(Si, Sj) to be the probability of a g-spin ((Sj, Sj) on a 
cavity graph marginalising over all other g-spins in the cavity graph. 

To determine the partition function on the cavity tree we decompose the problem 
iteratively; the partition function on the factor graph J-j-^ can be written in terms 
of partition functions on descendant factor graphs {Fk^jlk G nn(j) \ i}; the notation 
means k are the factors neighbouring j, excluding i. Referring to the schematic of figure 
[3] and references [29], [30] , the recursion is 



P 



,^i(Sj,Si)(X Yl 

kdnn(j)\i 



*y ] Pk^j(Sk, Sj 



exp{-/3Hj(Sj,S nn(j) )} , (9) 



and 



£ n 

Si,Sj k£nn(j)\i 



^ ] Zk^jPk^j(Ski Sj) 



exp {-(3Hj(Sj, S nn(j) )} .(10) 



The multi-body interaction is kept in a general form, whereas the intersection constraints 
are implicit in the appearance of the same spin Sj in each g-spin probability. 

The three degrees of freedom describing the probability distribution (Q can be 
expressed generally in terms of a coupling and two fields 



h' h b 



(Sj, Si 



with 



P.j,hSM( S v S i) K ex P {P {JSiSj + hfs i + h h Sj) } , 



(12) 



and (3 introduced as a known scaling. The superscripts are used to distinguish the 
ancestor (forward) spin field h* and descendant (backward) spin field h b . 

The probability distribution recursion (Q parameterized by ( JTT1) may be 
summarized as a non-linear mapping between the three-component cavity fields denoted 
by 



h i 



T ■ h f h b 

J— >n j- 



(13) 



allowing fl9]) to be written 

hj^i = f({h k ^j\k G nn(j) \ i}, 1/P) . (14) 

Note that the graph dependence is only in the arguments of T, and not in the mapping 
T itself. We have written here explicitly the dependence on the temperature 1/(3, for 
the convenience of taking the zero temperature limit in the cavity method. 
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3.0.1. Free energy Using (jSJ), the free energy can be calculated from the partition 
function at vertex as shown figure [2(a) 

— * — * 

- (3Fij = log Zi^j + log Zj^i - (3F E (hi^j, hj-n) , (15) 

where we call the edge free energy 

-PFBfaX) = log (j2 Si>Sj P^S,, S^S,, Sj)) . (16) 

An alternative calculation of the partition function shown figure 12(b) is achieved by an 
expansion about factor i, first creating the factor graph Ti 



(17) 



z{t) = zit) =j2 n 

Si j£nn(i) 

From this we have an alternative definition of the free energy 

-m= l0 SZj-,i- PFvdK^j enn(i)},l/0,Hi) , (18) 

defining the vertex free energy 

-f5F v ih 1 ...h c , 1/P, H) = log^n [P Kc (S c , S )] exp {-/3H(S , {S u . . . , S c })}j .(19) 

The trace Tr^ is a sum over all spin variables, and H is the Hamiltonian function. 

In the tree the number of factors is N and the number of variables is M = N — 1 . 
By summing the Fi and subtracting Fij, we have an expression dependent only on h for 
every edge 



(N - M)F = J2 F v e nn(i)}, 1//3, -J^ F e (k-tf, ^ 



(20) 



Another informative quantity is the free energy change when merging Co- 



descendant trees to generate a new tree (jSJ). This can be written as a function of 

(21) 



the Cj — 1 descendant fields 



AF^ = T F ({h k ^\k e nn(j) \ i}, 1//3) , 
where we define the site independent function 

1 

fjriih . . . ho-u V/3) = logTr 5 J] [P Kc (S c , S )} exp {-/3H(S , {S h S c })} .(22) 



c=l 



3.0.2. Marginals and order parameters From the cavity fields {h} we can derive 
relevant marginals. Taking the fields convergent on a particular neighbourhood we 
can define the joint probability distribution over the variables about a vertex 



P{Si, S nn{i) ) oc exp <j /3 \ ^2 (Jj^iSiSj + h^Sj + h^Si) - H^Si, S, 



nn(i) / 



(23) 
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and a similar disribution may be defined for edges 

P(Si, Sj) oc exp {/3 ((J^ + J^SiSj + (hUj + h)^)Sj + (hUj + } • (24) 

A marginal probability for variable z, P(Si), can be found marginalising over either (1231) 
or ( l24l) . which are constrained to give the same result. Summing the marginals, we can 
find the order parameters and standard extensive observables. The magnetization is 



N 

1 ■ 

m F 



1 N 

vEEw)' ( 25 ) 



N 

i=l Si 

and we take the spin glass order parameter to be 

2 



QF = Yl ( E S i P ( S i) ~ m P ) 

< \s, / 



(26) 



Finally it is useful to define the probability that a particular edge is in a dimer state, 
represented by a dimer magnetization 

m ° = if E E w(s> ^) , (27) 

Let the free energy density be / = F/N, then the energy density is 
1 N 

e=r jy E E Hi(ShS n n(i))P{Si,S nn (i)), (28) 

and entropy density 

s = /3(e-/). (29) 

5.1. Spin- symmetric solutions and the two-state model 

The recursive decomposition on a tree (114"]) terminates in a set of single g-spin cavity 
graphs on the leaves. With no external field, a cavity graph consisting of a single g-spin 
is described by a symmetric field 

h k -+j = (-A, 0, Of ; Z k ^j = 4 cosh(A) . (30) 

The zero values, h* = h b = 0, reflect the spin symmetry of the Hamiltonian, and 
the mapping f fl4|) does not break this symmetry. The recursion (jl4"]) is then non- 
zero in only one component Jk-+j- In a symmetric solution, which we later also call 
a paramagnetic solution, whether a g-spin (Si,Sj) is dimer (Si = Sj) or non-dimer 
(Si 7^ Sj) determines the recursion properties; the sign of constituent spins is irrelevant 
to the thermodynamics. 

For a general graph we find every symmetric solution of the "four-state" model (the 
model) in one to one correspondence with a solution of a "two-state" model defined with 
two state g-spin variables {SV,- = ±1 : ij G E}, and by a Hamiltonian 



#({*}) = ^e( E + A E<v 



(31) 
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The topology of the two factor graphs is unchanged other than in the presence in the 
two state model of factors corresponding to local fields on the variables SV,-. The "hard" 
intersection constraints are notably absent in the two state model. The quantity 

Pj-ti(Sij) oc exp [pJj^iSijj (32) 

becomes the relevant recursively defined object in the Bethe-Peierls or cavity method. 
The connection between solutions of the two-state model and the four-state model is 

Pj-+i(Sj, SilJ 7 , 4 states) = ^Pj^i{Sij = SiSjiJ 7 , 2 states) . (33) 

Since the marginal properties, and consequently many extensive properties such as 
energy, can be written in terms of these cavity probabilities. We can observe that 
if we identify Sij with SiSj then local and thermodynamic properties derived from the 
marginals are equivalent. Furthermore we find the difference between the free energies 
is accounted for entirely by a constant entropic term (C — 2)/ (2/3). In this sense the 
properties of the symmetric solution, other than its susceptibility to symmetry breaking, 
can be understood in the context of the simpler two-state model. 

This equivalence is rather surprising given that the models are not identical. One 
can make a two-to-one transformation of variables from the spins Si to a set of dimers 
{SiSj}, but the new state space is subject to constraints such that 1 = YlifijeL^iji f° r 
every loop L. This set of linear constraints on the variables are not independent, but 
are reducible to a set of M — N + 1 ~ (C/2 — 1)N independent linear constraints, that 
depend on the details of the graph structure. It seems however, that these constraints 
are irrelevant to the moments of the symmetric solution, and affect the entropy only as 
a simple function of /3 (independent of graph structure). For now we say only that since 
the loops are long (by the locally tree like assumption), and correlations decay rapidly 
in pure states (the fixed point purports to describe a single pure state), the effect of 
these constraints may be pushed to the boundary with respect to any neighbourhood. 
Since the nature of the constraints is not to favour either dimers or non-dimers per-se, 
we may anticipate a self- averaging effect at the boundary, and as such no bias would be 
present at, or persist from, the boundary. As such it is not surprising they have a weak 
effect on marginals even at low temperature. This argument relies on the assumption of 
weak correlation between a neighbourhood and its boundary at distance log(iV), which 
is expected to be reasonable in a pure state (at a stable fixed point for the probabilities, 
stable also against symmetry breaking). The effect on the entropy of the change of 
representation is discussed further in section 16.41 

3.1.1. Boundary conditions and symmetry breaking For a given boundary condition, 
a well studied property on trees, related to stability analyses in random graphs, is 
the stability in the core of a large tree to perturbations on the leaves. Typically one 
can define a variable at uniform distance from all leaves, and judge the stability of its 
marginal towards perturbations of the leaves - this can indirectly establish the presence 
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of long range correlations, susceptibility and other useful thermodynamic properties. 
Infinitesimal perturbations on a field 5h T = (5 J, 5h\ 5h b ), allow a linearized description 

Shk^j • (34) 

To establish asymptotic properties on trees we are interested in whether these 
perturbations introduced on the leaves decay under recursion about a given solution, 
such as the symmetric one. The two-state model stability analysis can differ from the 
stability analysis of the symmetric solution of the full model in that symmetry breaking 
perturbations of 5h* , 5h b can be considered in the latter. 



kGnnt 



dT 



dh 



3.2. The zero temperature recursion 

We wish to study the limit of zero temperature, where some simplifications are involved. 
We assume a decomposition of the field into energetic and entropic parts that are 
assumed to evolve on separable scales in the limit of large /3, 

K = /^-i/? + 0(l//3 2 ) . (35) 
Separating the leading order term in temperature from (fT4|) gives 

h E - ±h s = f E ({h E }) - ^f s ({hf}, {hi}, l/p) , (36) 
and allows a simplified representation of the energetic recursion 

JU = J E ({h k ^}) = — A - ^ ^2 SiSjD(Sj, Si, {h k -+j}) (37) 



fc£S = W) = - \ E SMSj, Si, { W) (38) 

S{ , Sj 

h& = h b > E ({h k ^}) = h f k '% S i D ( S i> 5 <» { W) ' ( 39 ) 

where 

D(S , S c , {h, h c -i}) = s mm ^ J ~ ( ]T S k ) + ((A - J k )S - h\) S k \ . (40) 



v fc=l / k=l 



3.2.1. Restrictions on the space of energetic fields In Appendix B we show that for 
any graph the distribution of fields generated self-consistently by ( |37]) - (I3"9]) is bounded 
at leading order in The two-state model solution is described by Jj_^ (I3"2i) . and 

a simple bound on these elements is found given connectivity Cj and an unconstrained 
set of descendant fields, namely, 

Jj^i G [ min J E ({h k ^}), max J E ({h k ^})} = [-(Q - 1) - A, (C,- - 1) - A] . (41) 
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In most cases this bound can be improved recursively by constraining the descendant 
fields self-consistently. 

For rational values of A additional degeneracy of the energy levels can greatly 
modify the thermodynamic properties of the system. In the simplest case of integer A 
the energy levels are even integers on any cavity graph: self-consistency requires the 
energetic field components to be integers. Since only a small bounded range of energies 
is available locally, energetic fields on every cavity graph must be from within a finite 



set. As argued in Appendix B we require a set closed under recursion (l3~Tj)-(l39l) and 



inclusive of any boundary conditions. 
4. Special Ratios 

Consider the properties of a factor of connectivity C within a factor graph, and assume a 
ground state of zero energy. A global ground state is certainly achieved if at every factor 
a local ground state is achieved. To achieve zero energy about a factor ([3]) requires a 
specific number of g-spins to be dimers in the neighborhood, but does not depend on the 
sign of the spin enforced by the intersection constraint. In the zero energy solution this 
number reflects a balancing of the next nearest neighbor anti-ferromagnetic interaction 
with the dimerizing (A < 0, favoring Si = Sj) or anti-dimerizing (A > 0) tendency of 
nearest neighbor interactions. The number of dimers in the neighborhood of a factor 
required to achieve a local ground state may be written as 

N c ,x = argmin Xe{Qr .. iC} j (x - J • (42) 

The minimizing argument must be integer, and will be unique everywhere unless C — A is 
a positive odd integer. For a given A, the set of vertex connectivities that are unlocked, 
producing degenerate solutions to ( l42l are 

C(A) = {C|C-A = 2a; + l,xeNo}. (43) 

Some integer values of A, we call the special ratios, may produce this dimer- degeneracy 
of local ground states, labeled by two consecutive integer values of x. We will call a 
factor that meets the special ratio criteria an unlocked interaction, whereas any other 
interaction will be called locked |31j . 

We say that a ratio is special with respect to a graph G(V,E), if amongst the 
vertices V of the graph, a finite fraction (0, 1] are unlocked. The set of ratios that allow 
unlocked clauses in a graph is 

A(G) = {X\Ci = C(A) for some % G V} . (44) 

Therefore we can say that a special ratio graph, is one in which the coupling ratio 
corresponds to an element in the set A(G). By contrast a special ratio ensemble for 
large N would be one in which replacing G by a typical graph A is in the corresponding 
set. Amongst the graph ensembles we study, it is only the case of regular ensembles for 
which every vertex is unlocked, in every other case the fraction of locked and unlocked 
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will vary with mean connectivity. For clarity we can discuss the implication for a regular 
graph with all vertices of connectivity C. 

We present a mean field argument to indicate the significance of special ratios. 
Consider the effects of adding a spin i to the system. Assume that in the absence of 
spin % the neighbors have freedom to take either +1 or —1 independently to achieve a 
ground state: the total number of ground states before adding the spin is 2 . After 
adding a spin the number of ground states is reduced to N gs = J2 x =n c x C\/[(C — x)\x\] 
for each state of Si = ±1. The number of assignments consistent with a ground state 
changes on addition of the spin by a factor r gs = 2N gs /2 c . 




Graph connectivity, C Mean connectivity, C 

Figure 4. Number of constraints per degree of freedom in regular graphs of various 
connectivity, (left) Asymptotically the ratio oscillates between asymptotes for the 
cases of Nc.x determined by locked and unlocked interactions, (right) For small C 
and A the oscillating pattern is apparent in linear graph ensembles, ensembles with 
sufficiently small constraint to variable ratios are good candidates for extensive ground 
state entropy. 



The factor l/r gs measures the number of constraints (implied by zero energy) per 
degree of freedom, which is intuitively linked to the probability a zero energy solution 
exists. The factor is high in the following two cases as shown in figure HJ First, when 
C — 1 < |A| for most vertices, in which case r gs equals 2 l ~ c for most vertices. Second, 
when C is large r gs becomes 0(1/ VC). Especially relevant to the previous discussion is 
the difference between odd and even connectivity nodes for given (integer) A. In special 
ratio ensembles two terms contribute to the sum N gs , and the entropy of the system is 
relatively high. As will be later shown for linear ensembles, r gs > 1/2 appears to be 
required for a zero energy ground state of extensive entropy. 

This argument would indicate that, for integer A, the greatest flexibility of the zero 
energy state space follows a pattern determined by the connectivity composition of the 
nodes. For example, if |A| = 0,2 one can anticipate significantly higher entropy for 
a zero energy ground state in a graph with odd rather than even connectivity. The 
opposite is true for |A| = 1. This opens the unusual possibility that as A is varied the 
entropy may oscillate, as will the free energy at finite temperature. Similarly, in the 
linear ensemble the proportion of odd and even connectivity nodes varies cyclically with 
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C. One can speculate that as constraints are relaxed and enforced when A or C changes, 
a periodic pattern of phases may occur. Indeed we will show this to be the case. 

5. Small graph studies 

In this section we study the energy and entropy density of ground states for iV < 200 
numerically. The linear connectivity ensembles chosen for presentation have small 
integer valued ratios (A) and small mean connectivity (C). These features allow a 
confirmation of the trends identified in section HI Further results are presented in 
Appendix A 



5.1. Sampling methodology and algorithms 

We generate simple graphs of A" vertices and M edges by a variation on the configuration 
model [32] . The ratio of the number of edges to the number of vertices in every sample is 
restricted to exactly C/2, and graphs of size A^ < 40 are restricted to be connected (every 
vertex is reachable from every other vertex along some sequence of edges) . An exhaustive 
enumeration method was applied to determine all ground states for small graph samples 
(N < 40). Deterministic greedy search and stochastic sampling, using the extremal 
optimization (EO) method [33l E4], was applied to larger samples (N = 0(100)). EO 
is an example of a stochastic local search method, performing a biased random walk 
in the space of spin assignments to discover ground states. This method was chosen 
because of the similarity between our model and algorithmically challenging constraint 
satisfaction problems where it has been successfully applied. 

Exhaustive enumeration of ground states was implemented by a branch-and-bound 
algorithm, exactly determining all ground states. Our EO parameters are calibrated 
against these exact results for small systems, and by further self-consistent analysis. 
For final data collection we adapted an algorithm provided by Stefan Boettcher that 
had been applied to a sparse spin glass problem [31]. Three parameters r, T T , R T control 
the quality of the ground state estimate obtained and are given alongside the data 
in the corresponding figures. Parameter r controls the probability to take locally 
suboptimal search directions, in a manner intuitively similar to the inverse temperature 
in Monte Carlo sampling. Parameter i? T is the minimum number of independent searches 
per graph sample, and is adaptively increased where variation between search result 
outcomes is large (hard regimes). Each search is initialized at a point in state space 
sampled uniformly at random. Parameter T T controls the maximum number of spin 
reassignments T T (1000 + ((N + l)/5) 3 ) per search. The cubic scaling with A" is chosen 
self-consistently with the amount of time required to find the ground state, for near 
optimal choices of r, in hard instances. To adapt the algorithm to determine the entropy 
was found to require a sampling time scaling as A^ 3,5 for relatively small systems A" < 50, 
we attempted only ground state search for the larger systems. 

Throughout the parameter space for the system sizes presented we achieved good 
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results by EO for r < 2, a parameter range in which locally suboptimal search 
trajectories are explored. The parameter r = 1.6 provided consistently good results; 
T T or R T were chosen sufficiently large that only a very small fraction of ground states 
appeared to be incorrectly identified. As an extreme contrast we studied greedy search 
methods, including the case r — > oo (numerically approximated as r = 10) of EO. For 
t — > oo we descend the energy landscape until no downward steps are available. This 
works well if the landscape is smooth and has a unique minimum: more generally to 
compensate for the trapping of the dynamics we can consider small T T and large R T . 

In the regimes where the number of constraints per degree of freedom is smallest, 
and we achieve zero energy, search time is fast for all algorithms - even linear in system 
size. Where the ground state energy is non-zero the number of samples required to 
reach a ground state increases significantly for all algorithms tested, with greedy search 
failing dramatically. For EO, with a well chosen r, the search time is slowest precisely 
in marginal cases, at the threshold between zero and non-zero energy. 

5.2. Results 

Energy and entropy results for the sampling methods on graphs of size N = 32 and 
100 are presented in figure EJ As connectivity and A are varied a non-monotonic trend 
exists in energy and entropy, in agreement with the special-ratio argument of section 
HI The pattern becomes more pronounced as system size increases. For A = 1 we see 
the energy minimized, and entropy maximized, about even connectivity. For A = 2 the 
trend is similar but for odd connectivity, increasing to A = 3 the pattern again reverses. 
A well-known feature of frustrated Ising models is the possibility of extensive ground 
state entropy, and for finite systems we see parameterizations that allow large entropy, 
in correspondence with regimes for which the energy equals the known lower bound 
(zero) . 

For N = 32 we demonstrate the entropy for A = 1 ± 5X, with 5X chosen as 0.01, 
the dimer degeneracy in the ground states is broken, with a large reduction in entropy 
relative to the results for the special ratio A = 1. However, the energy level splitting 
at each factor is proportional to 5X, the energy curve follows closely the unperturbed 
result. 

For C — > 2 and |A| = 1 we have the energy approaching zero and the entropy 
density in excellent agreement with the asymptotic result for the 1-dimensional next 
nearest neighbor Ising model, where s = log(l + \/5)/2 [35]. By contrast for A = 2 
there is a difference in behaviour of the energy for N = 32 and N = 100 approaching 
C = 2 in figure [5j This is explained by a difference in the graph sampling method: for 
N = 32 we have rejected all graphs except those that are connected (we approach a 
result dominated by a single unfrustrated loop); for N = 100 we allow fragmented graphs 
(we approach a result dominated by disconnected frustrated loops). Asymptotically for 
C = 2 and either graph ensembles the energy density must approach zero, since the 
number of loops is at most 0(log N), and each will raise the ground state energy by at 
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Figure 5. Statistics over 100 graph samples of size N = 32, and 1000 graph samples 
of size N = 100, for ground state energy (upper) and entropy (lower). Error bars are 
omitted for the entropy where A = 1 ± 0.01 for clarity. 



most 0(1). In the limit C — > 2 additional excitations might be anticipated as a fraction 
of the number of sites of connectivity greater than 2. 

Also apparent in figure |5] are points where the energy becomes non-zero 
continuously, implying continuous phase transitions. In figure [6] we give the probability 
that a graph will have a zero energy solution. Indeed, the drop in the probability 
becomes increasingly abrupt with increasing system size about C ~ 2.69. Results are 
shown from the exhaustive method that visit all states, and the heuristic sampling 
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2.5 2.55 2.6 2.65 2.7 2.75 2.8 
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Figure 6. The probability that a sampled graph (A = 1) allows a zero energy solution. 
EO has parameters (r,T T ,R T ) = (1.6, 100,7), the greedy result is achieved by EO at 
parameters (t,T t ,R t ) = (10,10,300). With 100 graph samples for N < 40, 1000 
samples for N = 100, 500 samples for N = 200. 

method, EO and a greedy search method. To compensate for the fact that the greedy 
method will be frequently trapped we consider many random minimizations (small T T 
and large R T ), but results are poor. Close to the transitions, the failure of the greedy 
method, and increasing search times can be taken as evidence for the appearance of 
many local minima that trap the dynamics in the vicinity of suboptimal configurations. 
The observation that the phase transition phenomena is correctly captured by EO, 
whilst the number of states visited (including repeats) is far short of the full state space 
2 N , demonstrates its efficiency 

For 2 < C < 5, A = 1 and 2, we have found numerical evidence for transitions from 
satisfied (energy 0) to unsatisfied (positive energy) phases in every integer interval of C 
as shown in figures. Numerically we find the transitions to be sharply defined for system 
size N = O(100), as in figure EJ with one exception which is the transition in interval 
C G (2, 3] and A = 2. In this case there is strong variability between sample energies 
and no clear cut-off for the system sizes studied. An important difference between this 
transition and all others will become clearer in light of the forthcoming cavity method. 

6. Cavity method 

The Replica Symmetric (RS) solution is presented in the context of the cavity 
method [31 |36j. The replica symmetry broken framework is provided at the level of 
energetic considerations, in the zero temperature limit [U 128] . 
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6.1. The replica symmetric (RS) solution 

The cavity method exploits the fact that the next shortest path connecting any two 
neighboring g-spins, after breaking the immediate path, is infinite in the large system 
limit, and that instances where the cavity graph rooted in a g-spin are not locally tree- 
like are statistically negligible. Assuming a pure state describes the thermodynamics, 
the point to multi-point correlations decay rapidly with distance, and so the estimate of 
probability Pt-+j(Si, Sj) is independent of Pj^i(Sj, Si) and other locally disjoint cavity 
graphs. The recursion (fl4"|) remains valid with the factorization of probabilities implied 
by this asymptotic independence. 



6.1.1. RS recursion condition We can make use of the self-averaging property for large 
iV to determine the fields by only a local recursion [36] . In the recursion the fields which 
are recombined at a generic factor can be represented by samples independently from a 
common distribution P(h). Self-consistency of this distribution requires 

P(h) = i ( C j f[ [dh h P{h k )\ 5 (X - f({h k \k = 1, . . . C - 1})) \ . (45) 
\ fc=l I c 

In the case of a regular graph we can identify one solution as P(h) = 5(h — h*). 
By symmetry of the Hamiltonian we expect the fixed point h* to be a symmetric 
field. Finding an analytic solution for modest values of C and zero temperature is 
possible, more generally it is straightforward to determine a numerical solution for any 
C, temperature, and A. 

A general method to establish a RS solution is population dynamics [36]. This 
involves an iterative procedure over a set of M. fields defined by h POP = {hi, ■ ■ ■ , hj^}, 
after t iterative steps. The approximation to the probability distribution over fields is 

Pt{ ^ = ji E S(h-h'). (46) 

An estimate at step t + 1 is created by iteratively updating once every element in the 
set h POP according to (IT4"|) by a stochastic sampling process from within the set. After 
a short transient a stable solution is achieved - one in which the moments converge up 
to small statistical fluctuations, Ai being a suitably large number. Detailed issues of 



the initial conditions and convergence are dealt with in Appendix C 



6.1.2. RS free energy To construct the free energy we can calculate the additional free 
energy when creating a typical graph of size N + 1 from a graph of size N. The free 
energy shift averaged over the growth processes gives the free energy density estimate. 
This is done by averaging the additive vertex free energy term, subtracted by the average 
edge free energy term to account for double counting of the free energy associated with 
the edges. That is, as in the case of a tree, 

f = Y^P(C)f v (l/P,C)-^f E , (47) 

c 
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where 



c 



fv = ! f[ [dh k P{hk)] F v ({h k }, 1/0, H c ) , (48) 

J k=l 

f E = J dti 1 P(ti 1 )dh 2 P(h 2 )F E (ti 1 ,ti 2 ) . (49) 

This is analogous to fl20l , composed of vertex ffl9l and edge free energies f|T6|) ; H c is the 
Hamiltonian applicable to a factor of connectivity C. The population dynamics solution 
again requires a replacement of the integration processes in (HHl) and (T49]) by random 
sampling, which is subject to fluctuation effects related to the population size M.. 

6.2. RS solution stability 

A first stability test of a candidate solution P°°(h) at the steady state t — > oo, considers 
the evolution of a perturbation P°°(h) + e(P*(/i) — P°°(h)). The second distribution 
P*(/i) can be chosen to maintain the normalization but with the moments of P°°(h) 
distorted by small amounts. 

In our study we focus on inhomogeneous graphs. We require a method compatible 
with population dynamics to test the stability of the fields towards linear perturbations. 
We achieve this by adding small random perturbations to every field and consider their 
evolution with time through a joint distribution P*(/i, 5h), maintaining the marginal over 
the first argument as P°°(h). A special case of the stability analysis can be undertaken 



for regular graphs; this is discussed in |Appendix C.3 



In this scenario the joint distribution of perturbations and fields at time t + 1, 
P t+1 (h,5h), is found by updating the leading order terms ( H5l) . and the subleading 
order through Under iteration the joint distribution can be expected to converge 
towards a distribution concentrated on 5h = (0,0,0); otherwise if there is not a single 
pure state we expect the perturbations to diverge. This can be quantified by the ratio 
of the second moment of the perturbations after many iterations to the corresponding 
moment 

1 f dhdShPHih, 5h)) \(5h f ) 2 + (5h b ) 2 + (5 J) 2 } , . 

77 = lim -log^ =; — - — L . (50) 

/ dhd5hP°{h, Sh) [(8hf) 2 + {5h b ) 2 + (5 J) 2 ] 

We are interested in the large t limit, but numerically we take t to be 0(100). In 
population dynamics, we expand the set hpop to be composed of a field and perturbation 
pair, the second moment thereby being the statistic of samples 

hpop = {(h^hi) (h N ,5h N )} . (51) 



Special care should also be taken in the zero temperature population dynamics 
method where we have energetic fields 0(1), and entropic fields 0(l//3), and stability 
should be tested in both parts. However, a frequent manifestation that indicates the 
inviability of RS solutions is the divergence of entropic fields under recursion from 
random (symmetry broken) initial conditions. Non-convergence of the entropic part, 
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given that we expect continuity of the solution in is an indication that RSB is 
required. 

6.3. The one-step replica symmetry breaking (1RSB) solution 

The next level of approximation in the framework is the one-step replica symmetry 
breaking solution. Our goal is principally to understand the critical behaviour 
associated with the competition of nearest and next nearest neighbor interactions at 
zero temperature. For this reason, as well as for issues of computational complexity, we 
focus on integer values of A, and linear connectivity ensembles with C > 2. This section 
gives a sketch of the methodology which has been presented systematically many times 
for generic locally tree like problems [28] . A standard stability analysis is considered as 



part of Appendix B 



In the method we restrict attention to energetic considerations at zero temperature, 
ignoring terms of 0(1/ 0) both in the consideration of the distribution of pure states, and 
by extension in the fields. The energetic field can be assumed to have integer restricted 
components for each pure state [281 EZ] , as was the case for a tree. 

Under the 1RSB assumption we must consider that at every site there exist many 
pure states. For every cavity graph the field, sufficient to describe a single pure state, 
must be replaced by a probability distribution over pure states V(f). At the energetic 



level, as shown in Appendix B , a finite set T = {h a \a = 1 . . . |T|} is able to describe 
the space of all realizable fields reproducible under recursion (I36I) . and is sufficient to 
describe the support of any solution. For our recursion we are required to describe in 
principle a field for every pure state, and the significance of each pure state should be 
weighted by its free energy. 

The fundamental assumption in the energetic 1RSB method relates to the free 
energies of the pure states [28]. The set of lowest energy pure states are assumed to 
have energies independent and identically distributed according to a Poisson process, 
with density 

P(E a ) = exp{fi(E a -E ref )} , (52) 

E re f being a normalizing constant. Equivalently we can anticipate the number M(e) 
of pure states at a particular energy level e to be exponential in the system size, and 
described by an exponent that must be a convex and monotonic function of energy 
called configurational entropy or complexity 

E(e) = ilogjV(e). (53) 

This gives a functional interpretation for fi as a variational parameter within the free 
energy, introduced by a Legendre transform enforcing the identity fl53|) . The complexity 
can be interpreted as the entropy of locally stable states at a particular energy density. 

Under recursion we expect that the pure states are each reweighted by the free 
energy shift they induce self-consistently with the Poisson distribution. It is argued that 
since in a pure state the fields on descendants are at long distances and uncorrelated, the 
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energy shift will equal that of a tree (j2ip . Therefore we should correct every mapping 
by a term with exp{— fiAF}, where fi acts now as a temperature over the pure states. 
Although the insight lies at the level of pure states, it is convenient to represent the 
properties by a field distribution without loss of generality. Since we have a finite set 
of fields the distribution can be represented by a vector V. Each component V a can be 
interpreted as the probability of a given field as being a (weighted) average over pure 
states on the cavity graph, 

U1 



V(h) 



0=1 



VJ(h a - h) 



(54) 



A relatively simple 1RSB solution might be obtained for a regular graph. Since 
we are considering the large N limit, the distribution of pure states is dependent only 
on the local properties of the cavity graph. In the regular graph every cavity graph is 
locally identical (tree-like), and so one case of ( 1541) describes every neighborhood. A 
self-consistent relation to describe this probability distribution over pure states is then 



V = T, 



RSB 



{V 



Ml 



(55) 



where using the free energy shift Tp in the zero temperature limit /3 — » oo applicable to 
independent fields ( |22|) . the components of the new vector are determined by 



>(c) 



l[h a = f({h ai ,...,h ac _ 1 }) 



(56) 



TrsbA^}) oc r£=Y [Eo c Li v£ 

x exp ^-/if F ({h ai ,h ac _ 1 },0) 

The indicator function I evaluates to one if the condition is met and zero otherwise. In 
the regular graph we must have all the distributions on left and right identical. 

In inhomogeneous graphs we must take into account that different cavity graphs 
will result in distinct distributions over pure states. Hence, the parameter relevant in 
the cavity method becomes a normalized distribution over cavity field distributions. 
Numerically this function can be represented by a population of M. vectors Vpop, by 
analogy with the population of fields used within the replica symmetric approach 

P ^ = Ji E KV-V'). (57) 
v'eVpop 

At any site we can consider independent samples from this distribution to reflect pure 
state distributions incident upon any vertex in the graph. The recursion (156]) becomes 
a componentwise mapping, whereas we would aspire to undertake an integral, in the 
population dynamics method we evaluate the integral approximately by sampling 

. c-i 



p(p) = / n [dv^pcp^ 

J c=l 



6[V a -T RSBia ({P^}) 



(5* 



a=l 



The total free energy is again given by ( 1471) weighted according to graph ensemble 
parameters, but with the 1RSB expressions 

c c 



MM/0) = Urn /n[ rfP ^ C )] lo Sn E P «c exp{-/i/3F y ({K ac },l//3)},(59) 
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replacing (Pig]) and 




dP{V')dP{V*) log ^ P a \P a 2 2 exp {-tfF E (H ai , h a2 )} , 



(60) 



replacing (Pj9]). Again we are only interested in the limit 1//3 = which restricts the 
argument in the exponent to integer multiples of /i. The variational dependence on /3 
from which energy and entropy are determined through the term /y (//, 1//3), is effectively 
replaced by a variational dependence on \i in both parts. From this expression we can 
construct the energy and complexity [2"5] 



Thermo dynamically meaningful solutions for \i exist where £(//) is non-negative, and a 
convex function of e{p). 

6.4- Cavity method for locked constraint satisfaction problems and b-Matching 

Ground states of energy zero in our g-spin framework can be considered as solutions to a 
constraint satisfaction problem (CSP). In the CSP framework each additive term in the 
Hamiltonian (JHJ) is considered to be a constraint. Local assignments producing minimum 
energy solutions of ([3]) satisfy the corresponding constraint, whereas realization of any 
other excited energy level is a violation. Globally a g-spin assignment is satisfied only 
if the energy is zero. 

By solving the 2-state model, we obtain a spin symmetric (or paramagnetic) solution 
to the full model, and we can show the 2-state model is equivalent to a well studied 
class of constraint satisfaction problems known as b-Matching [38] [39]. 

In the b-Matching problem, the task is to determine a set of edge states such that 
at each vertex i exactly b out of the Cj edges incident on the vertex are in state 1, the 
rest being in state —1. By analogy each g-spin in the two-state model describes the 
state of edges as dimer or non-dimer, with the zero energy condition requiring b = Nc t \ 
dimers about each factor (vertex), with Nc,\ defined in (T42|) . b-Matching has a trivial 
solution for b = or C (|A| > C — 1), since then an optimal matching is simply all 
dimers or all non-dimers. Other cases are non-trivial and have only recently been solved 
on random regular graphs. Standard b-Matching applies where Nc,x contains a single 
element, therefore the results we discuss do not apply to special ratios (where degeneracy 
is allowed). 

For random regular graphs, it has been shown that the satisfied solution always 
exists asymptotically in N, implying that the ground state energy of the two-state model 
is always zero. This result is most easily derived as a simple caveat to the 'contiguity' 
phenomena in random regular graph theory |32j. Furthermore, the solution space in 
the b-Matching problem in random regular graphs is known to be replica symmetric at 
zero and finite temperature [38], [39] . At zero temperature, the fixed point for the cavity 
recursion is described in our notation by 



e(/i) = Q-fif ; 



(61) 




(62) 
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where x = Nc,\/C. The entropy density derived via the cavity method is then |3T] 

S = C~ l0g icx) + t xlog ^ + log(l - x)] . (63) 

It is a simple matter to map the fixed point (|6"2"|) through (I3"3l to obtain the four-state 
symmetric fixed point, and again calculate the RS entropy. The entropy density is found 
to be reduced by (C/2 — 1) log 2 and hence negative for all cases. 

The following argument provides intuition on the contraction of the solution space 
when transferring from the b-Matching problem to the Ising Hamiltonian problem at 
any temperature or in any graph ensemble, that leads always to negative entropy in 
locked problems at zero temperature. This transfer involves attaching an Ising spin to 
every vertex, such that dimers and non-dimers are each formed in one of two states. In 
the spirit of the cavity method, we consider the effect of adding a node to the graph. 
This increases the entropy per node by log 2. Since the solution must be compatible 
with the intersection constraints, we have to subtract the entropy due to the double 
counting of the entropy log(2) per link, or equivalently C/2 log 2 per node. Overall, in 
the four-state model space, the entropy per node decreases by {C — 2)/21og2, relative 
to the entropy of the two-state model space. 

The negative entropy indicates that the spin-symmetric solution cannot be the 
physical solution for any parameterization for which the corresponding 2-state model is 
a b-Matching problem on regular graphs. Whilst this result is stated only for random 
regular graphs, we find it to be quite generally true that the paramagnetic solutions are 
of negative entropy at low temperature if the interactions are locked, so that we must 
seek a spin-symmetry broken solution. 



7. Results of the cavity method 

We present briefly some replica symmetric solutions, which qualitatively reproduce the 
phenomena of the small system experiments. We discuss only the positive A and C > 2 
scenario, where only solutions of zero magnetization are found. 



7.1. The paramagnetic solution 

We will describe in detail the paramagnetic solution for the case of regular graph 
ensembles, the salient features being relevant also for the case of the linear graph 
ensemble. This paramagnetic solution for the case of regular graphs, by local 
homogeneity, is described by a single cavity field h, and is compatible with the solution 
of the two-state model (1321) . The solution is a fixed point of the cavity equation 
(114p . considering a single field we must solve a polynomial equation, the order of the 
polynomial grows linearly with C — 1. 

To determine solution stability we can use the more precise method outlined 
in Appendix C.3 again avoiding a population dynamics approach. The eigenvalues 



describing the evolution of the mapping perturbations provide information on the 



NNNI models on random graphs 



26 



stability towards spontaneous symmetry breaking (a linear instability) and towards 
replica symmetry breaking (a non-linear instability). Linear instability may arise for 
sufficiently negative A, where an RS ferromagnetic solution emerges continuously at low 
temperature. We present results only for positive A, where we find the solution is always 
linearly stable. The second type of instability is towards replica symmetry breaking and 
is relevant to positive A. 
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Figure 7. Solid lines are the critical curves for the paramagnetic solution for regular 
connectivity graphs C = 3, 4, 5, 6, curves emanate from A = C — 1 respectively. To the 
right of the curves the paramagnetic solution is unstable towards RSB, everywhere else 
the solution is stable. The paramagnetic solution is stable down to zero temperature 
for small A. Markers for C = 4,5,6 indicate the crtical temperature below which 
entropy is negative. For C > 5 there is a reentrant behaviour for A < C — 1, but 
limited to the negative entropy regime. Dashed lines indicate the asymptote (1651) . 



In figure [7] we present the critical curves determined by the local stability analysis. 
For comparison, we plot the critical curves that would be anticipated from models 
of equal strength in nearest neighbor interactions, but without next nearest neighbor 
interactions. For many standard random and sparse graph ensembles including the 
regular and linear cases, we anticipate that asymptotically in large A the critical 
temperature will approach from below the critical temperature of the Bethe spin glass, 
given by the solution of 

1 = f>(g) C %~ 1) tanh 2 (2/3|A|) , (64) 
c=i 

where P(C) is the connectivity distribution describing the sparse model. This is because 
when A is large, interactions between spins are dominated by nearest neighbor anti- 
ferromagnetic couplings. In random graphs the presence of loops with odd number of 
edges causes frustration, and spin glasses are formed at low temperatures [40J. There 
is however a perturbation on this result owing to next nearest interactions which is 
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calculated in Appendix D The critical temperature found by an expansion in 1/ A is 

2 / 3C - 4 \ 

Tc ~ atanh(l/ v / C r ^I) \ ~ 2^/C~^l) ' ^ 

in agreement at leading order 0(/3X) with f|64"|) and in good agreement with the curves 
of figure CD 

The curves approach A c = C— 1 at zero temperature. It is noteworthy that for C > 5 
the curve becomes multi- valued for some A, indicating a reentrant paramagnetic phase. 
However, as will be seen in figure El the entropy demonstrates that the paramagnetic 
solution at this part of the curve is not viable. It is also interesting to consider the 
stability as 1/(3 — > 0. We find that where A < A c either the solution is marginally stable 
in locked ensembles (the eigenvalue approaches 1 from below in the zero temperature 
limit), or stable in special ratio ensembles. For A > A c we find the solution is unstable, 
whereas for A = A c (which are also special-ratio ensembles) the solution is stable for 
C < 4 and unstable for higher connectivity. 



Entropy density, s 
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Figure 8. Continuous lines indicate the entropy along the critical curves of figure [3 
curves corresponding to graph connectivity C = 3, 4, 5, 6 emanate from A = C — 1. 



In figure [8] we present the entropy of the paramagnetic solution at the critical 
temperature (s) indicated by figure [7J Along the critical curves we can have either a 
negative or positive entropy solution. For high enough temperature the local instability 
applies to a positive entropy paramagnet. This is consistent with the behaviour of 
nearest neighbor models, for which a continuous transition to a full-RSB phase is 
found [40] . The curve for C = 3 is exceptional in that everywhere the entropy is 
positive. So this curve may be describing correctly a continuous transition even at low 
temperatures. Otherwise at low temperatures, we have an unphysical behaviour of the 
entropy and must invoke a different phase to explain behaviour, nteractions. 

Despite the fact that for small A the paramagnetic solution is locally stable in the 
limit of zero temperature, we find that only if A and C are small, and A is a special ratio, 
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can we have a positive entropy solution. In table [TJ we give an exhaustive list of those 
regular graphs with C > 3 for which entropy of the paramagnetic solution is positive at 
zero temperature to 3 significant figures 

Thus amongst the locally stable paramagnetic solutions for regular graphs, it is only 
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Table 1. Zero temperature entropy of the paramagnetic solution on regular graphs, 
all positive entropy cases. 



a small subset of special ratios for which the solution is viable. The subset identified 
is consistent with the zero energy and large entropy results found in the numerical 
studies of section |5] For the majority of ensembles in which the paramagnetic solution 
is locally stable, but the entropy is negative, we must expect an alternative solution as 
the temperature is lowered. 

For inhomogeneous graph ensembles, the paramagnetic solution and its local 
stability have also been studied by means of population dynamics. The features found 
for the regular case seem to hold more generally. We find for A < C— 1 the paramagnetic 
solution is locally stable down to very low, or effectively zero temperature. By contrast 
for A > C — 1 we have local instabilities. As for the regular graph ensemble, only a 
subset of special ratio ensembles have positive entropy locally stable solutions at low 
temperature. 

1.2. RS spin glass solutions at finite temperature by population dynamics 

The paramagnetic solution is the stable and unique solution for all parameterizations 
at sufficiently high temperature. Initializing population dynamics with symmetric fields 
we always converge to this solution. Another solution at finite temperature has broken 
spin symmetry, and we call it the spin glass solution since the magnetization is zero 
and the spin-glass order parameter (1261) is non-zero. As the temperature is lowered the 
paramagnetic solution may be locally unstable towards a spin glass solution, or both 
solutions may coexist at low temperature. 

The extent of the spin glass phase is shown in figure [9] for a range of parameters. 
Above the critical temperature, populations of iV = 10000 fields converge to the 
symmetric solution within O(100) iterations from random symmetry broken initial 
conditions. Very close to the transition longer timescales are relevant, so that the 
critical temperatures presented in the figure should be considered upper-bounds, but 
still correct quantitatively up to O(10~ 2 ). 

The most interesting phenomena appear in the range |A| < C . This is the regime 
where nearest and next nearest neighbor couplings have comparable strength and the 
behaviour is influenced by proximity to special ratios. In the case of integer A we observe 
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Figure 9. The critical temperature 1//3, below which a spin glass solution exists for 
the RS equations is shown as a contour map. Parameters are sampled on a grid with 
increments AA = 0.05, AC = 0.05, the critical temperature found is represented by 
a grayscale: 1//3 = (black) to 1/(3 — 1 (white), contours are fitted by linear splines 
at the levels indicated by the legend. For some combinations of integer A and C the 
spin glass solution is absent at all temperatures, we assign the critical temperature to 
zero in these cases. A periodic pattern of critical temperatures is apparent, as A + C 
approaches odd integer values we find local minima in the critical temperature. 



that a paramagnetic solution can be the unique solution even at zero temperature. Two 
cases are shown in figure [TUl For A = 1, the zero-temperature paramagnetic phase exists 
when even connectivity dominates, agreeing with the prediction of (JHj) and reproducing 
the periodic pattern that characterizes the ground state energy in figure |5j The scenario 
for A = 2 is similar, except that zero-temperature paramagnetism exists when odd 
connectivity dominates. 

Figure [TT] demonstrates the extent of the RS spin glass solution in regular graphs. 
For non-integer A, we have only locked interactions and find the spin glass phases exist 
below some finite critical temperature, this temperature being lowest around integer 
values. The critical temperature is locally minimized at special ratio ensembles, and 
maximized at approximately their midpoints, following an odd-even pattern. Due to 
the special ratios there is no spin glass solution even in the limit of zero temperature, 
for some combinations of integer A and small C. 

Not shown in figures |9~HTT1 are the cases of larger C and larger A, where the critical 
temperature is typically much larger. As for regular graphs, we find evidence that low 
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Figure 10. Spin glass phase properties for (a) A = 1 and (b) A = 2. Data is 
sampled with increments A(l//3) = 0.01, AC — 0.05 and contours are fitted by 
linear splines. The uppermost lines indicate the maximum temperature allowing a 
spin glass solution to be found from random initial conditions of population dynamics 
(an RS spin glass existence curve). Solid lines indicate a discontinuous transition in 
qp, dashed lines indicate a continuous transitions from and a corresponding local 
instability of the paramagnetic solution. The solid (red) lines indicate contours in qp 
(|2"61 ordered monotonically decreasing in temperature, qp — 0.96, 0.92, 0.88, 0.84, 0.8. 
Labeled beside the x-axis are thresholds obtained from a zero temperature analysis, 
discussed in the main text, in apparent agreement with the 1//3 — > limit of the spin 
glass existence curve. 
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Figure 11. Spin glass phase properties for regular graphs C 
Data is sampled with increments A(l//3) = 0.01, AA = 0.05. 
as figure [TUl 



= 3 (a) and C = 4 (b). 
Lines have same sense 
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or zero temperature paramagnetic solutions are inviable either due to local instability 
(C > A — 1), or unphysical values for quantities such as entropy For A large relative 
to C we see again the emergence of a continuous transition following a critical curve 
consistent with a model dominated by nearest neighbor couplings as in ([Mil . 

It is noteworthy to distinguish between continuous and discontinuous spin glass to 
paramagnetic transitions within the figures, the spin glass order parameter either drops 
to zero discontinuously from a value qp ~ [0.7, 1], or goes continuously to zero as the 
temperature is increased. In the vicinity of discontinuous transitions a coexistence region 
exists between the paramagnetic and spin glass solutions. We observe that continuous 
transitions can occur for either: C ~ 2, i.e. models dominated by chains; or A •> C — 1, 
i.e. models dominated by nearest neighbor couplings. Elsewhere the transitions were 
found to be discontinuous. 

The RS population dynamics method for inhomogeneous graphs shows a similar 
local instability of the paramagnetic solution at large A to that uncovered in the analysis 
of regular graphs. In addition it demonstrates the discontinuous emergence of a replica- 
symmetric SG solution as the temperature is lowered or the ratio of locked to unlocked 
clauses is changed. The phase boundaries match qualitatively the mean field arguments 
of section H] and the numerical results of section [51 However, the RS spin glass solution 
is found to be locally unstable everywhere towards further levels of replica symmetry 
breaking, the specific heat and entropy behave unphysically in all spin glass solutions. 
Furthermore, even where the paramagnetic solution is locally stable, it is frequently 
found to be a negative entropy solution below some critical temperature. Higher level 
RSB is required for a correct description; but we develop the cavity method henceforth 
only for the limit of zero temperature. 

7.3. Replica symmetric solutions at zero temperature 

In the zero temperature limit the RS population dynamics involve energetic 0(1) and 
entropic 0(l/f3) fields. If these fields are initiated symmetrically, the population can 
converge to a paramagnetic solution. The solution can be locally stable, and of positive 
entropy close to the special ratio values. However, as for finite temperature, the 
paramagnetic solution is often of negative entropy away from these points. 

We can by contrast find a critical behaviour if the fields are initialized as 0(1) 
random numbers. The solutions are then paramagnetic ones only for odd (even) integer 
A, and with sufficiently many even (or odd) connectivity vertices. We mark in figure [TUl 
the critical values in C at which the paramagnetic solution becomes the unique solution 
achieved from all initial conditions. These points are in good agreement with the limits 
of the finite temperature results. 

Outside this regime where the paramagnet appears to be the unique solution, we 
find the entropic fields diverge under iteration so that no solution can be found, with one 
exception. An RS spin glass solution is found for special ratio A = 2 and C G (2.4, 2.49]. 
This solution has spin-symmetric energetic fields, but spin-symmetry broken entropic 
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fields; qp is non-zero, and goes continuously to zero at C = 2.49. Although population 
dynamics converges in this case, this RS spin glass solution is locally unstable to RSB, 
like the finite temperature spin glasses. 

7-4- Energetic 1RSB solutions 

To gain insight beyond RS we look to 1RSB and focus on graphs with low connectivity 
and integer A. We will skip the discussions on the cases of non- integer A, large integer 
A, or large C, since the state spaces in those cases are strongly constrained, as argued in 
section H] and so far we have not found evidence of zero temperature phase transitions 
in these regimes. 

To facilitate the calculation of the energy and complexity in ( I6ip . we first compare 
the behaviour of the full Hamiltonian ([3]) and the Auxiliary Hamiltonian ([6]). Results for 
the complexity E(e) are shown in figure [121 The new model provides a lower bound for 
the energies at given complexity, and an upper bound for the complexity at given energy, 
where it exists. There seems to be a good overlap in these two parameters so that it is 
useful. However, the absolute values predicted for the ground state energy, and spinodal 
point (solution of maximum energy), differ quantitatively and to consider a better 
auxiliary, for example we may approximate the local energies by the set {0, 1, oo}, where 
infty corresponds to forbidden energies when the local energies of the full Hamiltonian 
are 2 or above. As opposed to {0, 1}, this might tighten the curves for such a purpose, 
at the cost of complexity. We have also compared the results for non-integer C, the 
agreement between the two Hamiltonians becomes more convincing. This is because 
figure [12] presents the case of maximum frustration, in which local excitations (where 
the Hamiltonians differ) are necessarily realized. In regimes where the ground state 
energy is close to zero the gap between the curves becomes correspondingly tight. 

The complexity curves in figure [12] are typical of the regimes of locked constraints, 
say (A, C) = (1,3) or (2,4). In these cases E(e) starts from a negative value at e = 0, 
indicating that the zero-energy state is thermodynamically negligible. At this point the 
configuration parameter fi approaches infinity. When e increases, S(e) increases and 
the corresponding value of ft decreases from infinity. The intersection with the energy 
axis yields the ground state energy at zero temperature. When e increases further the 
positive values of S(e) ends in a cusp, indicating the absence of metastable states of 
higher energy. This point indicates a dynamical transition, the spinodal point. This 
means that for a dynamical process in which the temperature is lowered the state will be 
trapped in metastable states at this energy due to their configurational entropy. Figure 
[12]also shows the upper concave segment of the complexity curves, which are unphysical. 



The stability analysis of 1RSB for regular graphs is discussed in Appendix C and 
follows standard local tests called type I and type II [3]. The solutions described by the 
curves in figure [121 are found to be type I stable throughout. Testing type II instability 
it is found that they are unstable everywhere: both the metastable states, and the 
thermodynamic state, must be considered within a higher level of RSB to correctly 
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Figure 12. Complexity curves for the locked problems on regular graphs: C = 3, A = 1 
(top) and C — 4, A = 2 (bottom). Physically relevant solution for complexity is the 
convex upward part from the point of maximal complexity to the intersection at zero 
complexity. Phenomenology of the full © and auxiliary © Hamiltonian solutions are 
similar, but with significant errors in absolute terms. 
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Figure 13. Zero temperature complexity curves under the assumption of zero energy, 
A = 1 (top) and A = 2 (bottom), and a linear connectivity ensemble. Dynamical and 
thermodynamic transitions values are indicated by labeled arrows. 



describe all system properties. 

To locate the transition points between spin glass and paramagnetic solutions, we 
plot E(0) as a function of C in figure [TBI Since energy is restricted to zero, the excitation 
structure that distinguishes the auxiliary and quadratic Hamiltonians are irrelevant, 
the results are the same. Since E(e) is an increasing function, negative values of S(0) 
indicates a higher ground state energy e is relevant, and positive values of E(0) indicate 
the presence not only of a zero energy ground state solution, but of exponentially many 
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metastable spin glass states in the system. The intersection of E(e) with the C axis 
indicates the transition point between spin glass and paramagnetic phases, and the 
points where the positive complexity disappears indicate dynamical transitions [28j H] • 
These transitions are labeled in figure [T3j 

A sequence of transitions can be observed when the value of C alternates between 
odd and even dominated regimes. The curves of S(0) have minima in the regular graphs 
with locked interactions, such as the point (A, C) = (1, 3) indicating that the zero energy 
states are absent. When one moves away from these points, the inclusion of unlocked 
interactions into the graph causes S(0) to increase until it becomes positive at the 
transition points. Comparing the RS and 1RSB predictions, we observe that the 1RSB 
results yield a broader spin glass regime. For example, around (A, C) = (1,3), the RS 
and 1RSB spin glass regimes are given by 2.76 < C < 3.33 and 2.69 < C < 3.45 
in figures HD1 (a) and [13] (left) respectively. Besides the thermodynamically stable 
phases, the 1RSB analysis also predicts the existence of metastable spin glass phases 
at 2.65 < C < 2.69 and 3.45 < C < 3.56. Comparing with the numerical results in 
figure [6] we can see good agreement between the thermodynamic transition prediction 
and the point at which P(e = 0) goes to zero with increasing system size; metastable 
states would also explain the algorithmic slowdown of stochastic local search methods 
approaching this transition. 

When C increases we observe that the paramagnetic phases narrow down until 
the entire paramagnetic phase coexists with a metastable spin-glass phase and there is 
no longer any dynamical transition. An example can be found at 4.83 < C < 5 for 
A = 2. When C is sufficiently large, we expect that the paramagnetic phase at zero 
temperature disappears due to the increasing number of constraints, and E(0) becomes 
negative everywhere. Also for C < A we find that E(0) increases continuously from 
negative to zero without a dynamical transition. This is apparent at (A, C) = (2, 2.4) 
in figure [13j We have only shown results for A = 1 and A = 2, similar properties may 
characterize other positive A. 

In figure HH with A = 1 we find the ground state energy for the 1RSB solution, given 
by the energy at which the complexity first becomes positive on increasing /i. Some 
granularity in \i of size 0.01 is visible in more computationally challenging regimes, 
leading to the appearance of step-like artefacts. Again the difference between the 
auxiliary Hamiltonian and full model is small. The curve for small C indicates a very 
low value for the energy by contrast with other curves, only 0(1/100) factors are excited. 
A second spinodal line, tracing the structure highlighted by figure is not drawn but 
follows a similar trend. Comparing the energy to early experimental results, there is a 
systematic trend apparent in figure [5] that the energy is decreasing with system size. As 
well as being qualitatively similar to the numerical results, the curves in figure [141 seems 
to be a lower bound to the finite system experiments as anticipated. 

The analysis of the 1RSB and RS energetic fields, at zero temperature, indicate a 
spin glass transition at A = 2 and C = 12/5. In the full RS analysis, considering also the 
entropic contributions, a spin glass transition at C = 2.49 is observed. We argue that 
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Figure 14. Ground state energy for the linear ensemble with auxiliary and quadratic 
Hamiltonians (C < 3). (Top) A = 1. (Inset) Close up view of the region near (7 = 3. 
The configurational parameter relevant approaching C — 3 is marked for the two 
Hamiltonians. (Bottom) A = 2, (Inset) Close up view of the region near C = 2.4. The 
transition at 2.4 is due to a percolation transition and data is noisy, the difference 
between auxiliary and quadratic Hamiltonians is small in this regime. Data is not 
collected in 2 < C < 2.03 for numerical reasons, the curve is anticipated to continue 
smoothly to (C,e) = (2,0). 



the value C = 12/5 has a topological significance and probably represents the correct 
transition point from a zero to non-zero energy solution. This value corresponds to a 
percolation threshold. For C > 12/5 there begin to exist percolating clusters in the 
factor graph on which all factors have connectivity 3. Symmetric fields are generated 
and propagated by such a substructure. Conversely, for C < 12/5 the graph is composed 
of connectivity 3 factors interspersed by many factors of connectivity 2. The interactions 
along these chains restrict g-spins to non-dimer states, then for a given configuration 
of g-spins on the boundary, all internal spins become aligned anti-ferromagnetically. 
In a random topology where this substructure percolates, many closed paths lead to 
frustration, a raised energy and spin-glass behaviour. The details of this argument help 
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to explain the exceptional regime of strong finite size effects seen in experiments (section 

ED. 

8. Discussion 

We have studied the Ising model with nearest and next nearest neighbor 
antiferromagnetic interactions on random graphs with varying dominance by odd or 
even connectivities. The model is closely related to many problems in optimization, such 
as maximizing the diversity of colours on random graphs and matching. Experimental 
observations are made by exhaustive branch-and-bound search for small systems and 
EO for larger systems. 

Theoretical analyses were undertaken by the cavity method in both the RS and 
1RSB approximations. The results demonstrate the existence of both paramagnetic 
and spin glass phases at zero temperature. Zero temperature paramagnetism can be 
explained by the high configurational freedom in special ratio ensembles due to unlocked 
interactions, and the transition to spin glasses, at small A and C, can be attributed to 
the inclusion of locked interactions as the average connectivity moves away from the 
special ratios. 

To deal with nearest and next nearest neighboring interactions on equal footing, we 
have developed the formalism of generalized variables. This enables us to conduct the 
thermodynamic analysis based on locally tree-like recursions, resembling those used to 
study the equilibrium and dynamical properties in the Bethe lattice and Cayley tree. 

It is also worthwhile to consider extending our algorithmic strategy to other 
computationally hard problems. We first perform simulations on small systems using 
exhaustive enumeration, which already reveals the effects of special ratios. Since the 
system sizes they accessed are too small to accurately estimate finite size effects, we use 
in larger graphs the heuristic EO. The transition point found by EO is increasingly well 
defined with increasing system size. Its power is evident in its ability to discover ground 
states for system sizes up to 0(100) in a short computation time. 

As far as we know, this study is the first in applying the 1RSB cavity method to Ising 
models on random graphs with next nearest neighboring interactions and inhomogeneous 
connectivity. However, both the RS and 1RSB solutions are not sufficient to exactly 
describe the spin glass phases at low temperature. The description of the energetic 1RSB 
approximation, which has been studied at high numerical precision in many graphical 
models, becomes quite complicated owing to the multiple components of the cavity fields 
(3 for h in our case). Nevertheless, the approximations are in qualitative agreement with 
the simulation results for small graphs. 

Looking forward to studying models with higher connectivity, greater distance 
interactions, or Potts spins, it will be necessary to introduce simplified approaches 
to keep the problem tractable. Our work also demonstrates some useful directions 
along this line. Applying the population dynamics in the 1RSB framework we have 
demonstrated an efficient method to enumerate the support, in the space of fields, for the 
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order parameter. The use of the auxiliary Hamiltonian is also able to produce accurate 
results on the transition points, and reasonable estimates at finite temperature and 
finite energy. Other simplification strategies, not yet explored in this paper, include the 
Gaussian approximation for high connectivity graphs and message passing algorithms 
approximated with reduced number of variables [21]. In exploring the possibility of 
Potts spins, a proposal simplifying the complexity of calculating the marginals has been 
recently proposed, using a belief propagation algorithm on generalized states [25] . 

We have studied regular graphs and inhomogeneous graphs with a combination of 
connectivity \_C\ and \C] . We expect that in other ensembles the prevalence of odd 
or even connectivity nodes will remain a dominating factor in determining the phase 
diagram, provided that the connectivity is not too large and A is small by comparison 
with vertex connectivity. For example, if A = 2, zero temperature paramagnetism will 
only be found in graphs with a significant fraction of odd connectivity nodes. However, 
in the cut-Poisson ensemble, which has a nearly uniform distribution of odd and even 
connectivity nodes for any C, we do not see the periodic trends in free energy observed 
in the linear ensembles with variation of C and A. We anticipate a similar pattern for 
other standard connectivity distributions. 

Special ratios should also affect the low temperature properties of graphs with A 
close to integer values. Since A is not an integer, genuine unlocked interactions no longer 
exist. However, since A is close to integer values the energy gap to the first excited state 
is small. Thus the special ratio effects are manifested in a reduced critical temperature. 

Besides the Ising models, special ratio effects are also expected in Q-state Potts 
models. For example, in the colour diversity problem, which corresponds to A = 1, 
the local configuration for Q = 4 of a node is locked at C = 3, and unlocked for 
C = 4, giving rise to a spin glass to paramagnetic phase transition when C increases 
from 3 to 4 in a linear graph ensemble [23J. In general, for A = 1, the multi-body 
interaction is locked when C + 1 is an integer multiple of Q, and the corresponding 
ground state local configuration on g-spins translates to each of the Q colours appearing 
(C + l)/Q times, with a total of C + 1 colours decorating a node and its nearest 
neighbors. Indeed, we have calculated the ground state entropy in the RS ansatz for 
A = 1 and 2, and found that they are particularly high in special ratio ensembles. One 
may wonder whether transitions between zero temperature paramagnetism and spin- 
glass phases can be observed with variation of A or C for some ensembles. However 
for larger Q, the special ratio ensembles are restricted to relatively high connectivity, 
with many constraints per degree of freedom. Even for those special ratio ensembles of 
lowest connectivity the paramagnetic solution becomes unfeasible, the presence of next 
nearest neighbor interactions imposes stringent constraints on the colour configurations 
that renders the entropy negative. Nevertheless, the special ratio effects may again be 
exhibited in a modulation of some thermodynamic quantities with variation of A or C, 
such as the dependence of spin glass transition temperatures on the average connectivity. 
This is a challenging topic for future studies. 
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Figure Al. (a) Mean entropy density (s) and (b) mean energy density (e) for linear 
ensembles. Results of (s) and (e) for Cut-Poisson ensembles are shown in (c) and (d) 
respectively, 100 graph samples in every case. Error bars are typically smaller than 
the point size. 



Results are presented to complement those discussed in sectional Figures lATl fa) and 
(b) demonstrate the entropy and energy statistics for linear ensembles of size N = 24. 
For A = 0, 3 the energy is minimized at the special ratios, while the entropy is maximized. 
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With increasing A, regimes of ground state energy zero become rarer, and unattainable 
above some threshold in C in line with the mean field argument, RS and 1RSB results. 

Figures lAlT c) and (d) show results for Cut-Poisson ensembles with N = 32. These 
figures demonstrate a nearly monotonic trend in entropy and energy statistics; strong 
mixing of locked and unlocked factors for all C prevents the appearance of modulated 
behaviour due to dominance by even or odd connectivity. 

Appendix B. Zero temperature field support 

Whereas the space of energetic fields is sampled statistically within the RS approach 
of the main text, our implementation of the 1RSB cavity method requires an exact 
description of the space. We establish properties of the space of energetic fields in this 
Appendix. Attention is focused on integer A, and we consider ensembles of minimum 
connectivity two so that boundary energetic fields can be ignored. This allows us to 
restrict attention to spaces of fields defined only by the zero temperature mappings 

dHZD-P}. 

The exact knowledge of the support allows a higher fidelity representation of the 
probability distributions manipulated in the cavity method for either RS or 1RSB. 
However, the methods outlined prove to be infeasible in many graph ensembles due to 
the restriction of studying ensembles with some maximum connectivity C max , and the 
unfavorable scaling of the space of fields with the connectivity. In principle a population 
dynamics approach can be applied to the 1RSB analysis for these cases [36J . 

Appendix B.l. Bounds on the space of energetic fields 

The components of the energetic fields, assuming that they are generated by recursion of 
the mappings (!371) - (l39|) from some initial condition, will be shown to be bounded. This 
reflects the intuition that changing the state of the g-spin at the root of the cavity tree 
allows only a restricted number of energy shifts in the ground state energy regardless of 
boundary conditions. 

For analyzing the zero temperature field support, it is more convenient to introduce 
a new representation of hf_^ to replace the representation ( Jj-w, rf^ii ^-h) usec ^ ^ n ^ ne 
main text. We can do this by decomposing the cavity probability according to 

Pj^i(Sj, Si) = Pj^>i(Si\Sj)Pj-ti(Sj) . (B.l) 

where Pj-n(Si\Sj) is the conditional probability of Si given Sj on the cavity graph 
J-j-H, and Pj^i(Sj) the marginal on the same graph. For Ising spins we can express 
Pj-+i(Sj, Si) in terms of three probabilities Pj^i(Si\Sj = +1), Pj^i(Si\Sj = —1) and 
Pj-+i(Sj). In turn, the fields and bj^i can be introduced to describe the three 
probabilities, namely, Pj^i(Sj) oc exp{(3bj_>iSj} and Pj^i(Si\Sj = ±1) oc exp{/3 fjL^Si} . 
Combining these cases, we can write Pj^ t [Sj, Si) oc exp[j3Ej-n], where in the limit of 
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low temperature 

iW%,ft) = - j£^ = l/jIJ ) 5, .(B.2) 

We henceforth work with fields only at the energetic level. Comparing with (11 II) we 
obtain the transformation in the limit of low temperature given by 

( l^j-ti + — — "Tj'-hI ) 

= frj-i + " 2 (B ' 3) 

= ^ ± J^i ■ (B.4) 

Appendix B.l.l. Bounds on fields describing symmetric solutions The symmetric 
solution is described with /r^ = — /t^, and = 0. At the energetic level 

the recursions become identical to those of a two-state model ( 133|) . the parameter 
Jj_u — ► ^fjt+i can be considered to describe either the two-state model or the symmetric 
solution of the NNN model. 2/3 becomes the logarithm of the ratio of likelihoods 
for dimer and non-dimer states of the ancestor g-spin given its descendants. 

A simplified form of ( 1371) applies for the two-state model, taking S = 1 or —1 to 
represent a dimer or non-dimer state 

Jj^t = J E {{J k ^}) = -A - ~ SD^ 2 \S, {J k ^}) , (B.5) 
where (l4"0l) becomes 



2 

s=±i 



)/C— 1 \ C— 1 C— 1 I 

^ \fc=i / k=i k=i J 

The first two quadratic terms arise from the next nearest neighbors interactions, and 
the final two terms are nearest neighbor penalties. A maximum and minimum value 
for an ancestor field Jj->i is found by freezing the descendant fields to uniformly large 
values of either all positive or all negative in sign. Doing this §k = sign(J fc ^) 

and maxima and minima of the expression are derived giving the bounds ( 14T|) . The 
bounds can be somewhat tightened by recursively reintroducing the bounds found on 
the ancestor, into the descendants. 



Appendix B.1.2. Bounds on fields describing symmetry broken states Without 
restriction to a symmetric solution we can present bounds on the components 
and £>_,•_►». We can take first the case of ff^i, determined as in (1B.2|) by 

/^ = -^^^(±l,^), (B.7) 

5, 
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where the energy is given by the low temperature limit of (I19p . 

Ej^Sj, Si) = min 5lv .. i 5 c _ 1 j (XSj + Ylk=i S k) S i ~ \ YSZx (fiUj + fk^j) S i 

+ Eti 1 [I (lAJUil + fLjSi - \f^j\ + + eETi 1 (b-8) 

+ ASj — bk^j] Sk} ■ 

From ( 1B.7I) we note that /-^ are determined by the energy change when the ancestor 
spin Si flips at a fixed value of Sj. Since from flB.81) the energy primarily depends on 
Si via the term (XSj + Yl Sk)Si> bounds for fp^ are determined in the cases that the 
fields . are strong, so that the values of Sk are fixed. Taking the combinatino of 
Sk that yields the maximum and minimum arguments, we find that /t^ and — /Z^ 
are subject to the same bounds as J^-h in (|4Tj) . For this representation of the cavity 
probabilities, /j^ are bounded the same way for symmetric (two-state) and symmetry 
broken (four-state) boundary conditions. 

By contrast is unbounded, given unbounded descendant fields. For brevity 
we demonstrate the case of uniform connectivity C on ancestor and descendant nodes. 
However, we can apply the bounds established for fjl^ to the descendant fields f^j- 
bj^i is given by 

b^i = ~\ Yl S i ™ in S^} . (B.9) 

Taking the limit that bj.-^.j are very strong so that Sk are aligned with bj^, we obtain 
for non-negative A 

c-i 

b j ^ i = Y,& j -K^-C+l)--[\2X-C + l-X\-\2X-C + l + X\] , (B.10) 

fe=i 

where X is the number of descendant nodes with Sk — 1. The upper bound is reached 
when X = and . = A + C — 1 for all k. This yields 

^ < (C-1)(C-1+A) + (C-1)A-^ [\C - 1 + A| - |C - 1 - A|] .(B.ll) 

The lower bound, as well as the case of negative A, can be derived analogously. 
Summarizing, 

\h \<Ti*(\\\r \- I (C-l)(C max + 2\X\)-\X\ if|A|<C-l 
\b^\<H(\X\,C max )-^ {c _ 1){Cmax _ 1 + m otherwise) 

generalizing from the case of uniform connectivity we have C max which is the largest 
connectivity in the graph ensemble, the extreme values of the bound are only achieved 
if all descendants are connectivity C max and the bounds (l4"Tj) saturated. 

A combination of these two bounds places a bound on the components of h, 
the domain is restricted to a polyhedron within any ensemble of graphs of maximum 
connectivity C max . The polyhedron has a simple description, but its volume can be 
compressed further by recursively reintroducing the bounds found on the ancestor into 
the descendants. In so doing the bounds involve a complicated coupling of {/ ± , b} in 
general, and are evaluated only numerically. 
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Appendix B.2. Integer A discrete set support 

Integer A cases are the focus of the main text. At these values the energy level 
differences are restricted to be integer. The field components, which represent these 
relative differences on cavity trees, are thereby restricted to integer values. Since the 
components are also bounded, only a short list of labeled fields can describe any possible 
pure state at any point in the graph, at the energetic level. Some small set of fields 
is sufficient to describe the support for the distributions describing either 1RSB or RS 
solutions. 

Our 1RSB method utilizes this discreteness, but involves summations over pure 
states. A greater variety of ensembles, and fidelity of results, can be obtained if small 
sets can be used. Our objective is to find the smallest set that describes correctly the 
solution support. 

A sufficient set, defined Yq a , is the set of all fields with integer components, 
consistent with the previously derived component-wise bounds. An analogous set is 
defined when restricting to the symmetric fields, ^c^f x- Both these sets can support 
solutions, and in some cases the same unique (paramagnetic) solution; however, it is 
much simpler to solve the recursions on the latter set, owing to its small size. The sets 
are completely defined by C max and A, their sizes are 

\ r clat\\ = (2H b (C max , |A|) + l)(2C max - l) 2 l r c°"" d A 2) | = (2C max - 1) . (B.13) 

Without loss of generality in the application of the cavity method, it is possible to 
restrict attention to self-consistent sets. A self-consistent set contains only fields that 
can be generated by some combination of fields in the set under mapping ()36l)-(l39i). 
The largest such set will certainly be no bigger than Y bound , and can be discovered by a 
recursive pruning process, removing irrelevant fields. 

Consider a mapping from a set V® to a set r( t+1 ) defined 

r(* +1 ) = T r (rW) = {h E \h E = f E ({hf}), for some {tif } G T®,C} , (B.14) 

where T E is the mapping defined by A and C (|36l)-(l39l). and C can be any connectivity 
within the graph ensemble. This mapping may be iterated and assuming convergence 
is achieved, r^ +1 ^ = the set derived will be a self-consistent one. 

Such an iterative procedure is not guaranteed to converge for an arbitrary choice 
of initial set, but taking = T b ^ und n ,, in effect all fields as the initial condition, the 
mapping is a contraction and must converge. For the linear ensemble with maximum 
connectivity C max , minimum connectivity C m { n and coupling ratio A, we obtain the set 
rl r , \ . For the symmetric set of fields we can similarly define T% n x, which 
is arrived at by the same process but taking r^ ^ = Y b ^ and ^ x y 

By construction T 4 is the largest set of fields self-consistent under the mapping Tp, 
it is generally a much smaller set than Y bound . The bounds of section Appendix B.1.2 



predict a scaling of 0(C 4 ) for the redundant set T R on which our most precise methods 
rely, and we are unable to study sufficiently large set sizes to establish a more favorable 
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scaling. By contrast T 2 has a simple description for the linear ensemble. It can be 
derived analytically, for ensembles in which C max < C min + 1, 

j . g J [~(^rnax — 1 — |A|),C moa; — 1 — |A|] if |A| < C m i n ~ 1 ^gx 
| [Cmin ~~ 1 — \X\, Cmax ~ 1 — if |A| > C m j n — 1 . 

The derivation of (IB. 151) is outlined as follows. The coupling is given by 

[Ej-ti(l, —1), Ej-yi(l, l)]/2. Let L and if be the lower and upper bounds of the couplings 
respectively. When the couplings reaches the bounds, it is reasonable to expect that 
the descendant couplings also reach the bounds. Suppose among the C — 1 descendants, 
m of them take the upper bound and n of them the lower bound. The cavity energy is 
reduced to 

E^(l, Si) = XS, + (Si + X)(x + y) + l -{x + yf - Hx - Ly , (B.16) 
where x = ^2 k \j =H Sk and y = J2k\j =l When A < C—l, empirical observations 

I k — yj I h — vj 

show that H > > L. Hence the minimum energy is given by 

( XSi + (H- L)m - ( L - x - s ^ iim^nKL-X-S, 
E ] ^ l {l,S l ) = \ XSi + (H- L)m - ( H - X 2 S ^ 2 if m-n>H-X-3i (B.17) 

[ XSi — Hm + Ln+(Si + X)(m — n) + ^-^- otherwise . 

Respectively, the above three cases are obtained at the x = m boundary, the y = n 
boundary and the corner (x,y) = (m,—n). These energy expressions result in the 
couplings summarized by 

—L if m — n < L — X 

Jj^i = { -X + m-n \iL-X<m-n<H-X (B.18) 
— H if H — X < m — n 

Since m — n is bounded by ±(C + 1) we see that H = —L = C — 1 — |A|. Generalizing 
the result to multiple connectivities, the first case of (IB. 151) is obtained. 

When |A| > C max — 1, L — H — for C = C max . This means that the next 
neighbor interactions (represented by C max — 1) are completely neutralized by the nearest 
neighbor interactions (represented by |A|). The energy minimum starts to be dominated 
by nearest neighbor interactions, with descendant spins aligning anti-parallel to Sj for 
positive A, and parallel for negative A. Hence we consider the solution with x + y = l — C 



for positive A and C — 1 for negative A and obtain the second case of (IB. 151) . 
The set size 

|r CiA | = 2m&x{C max - C min ,C max - 1 - |A|} + 1 , (B.19) 

increases linearly with C for fixed small |A|, and decreases linearly with |A| to 
2(C max — C m in) + 1- For large |A| the largest self-consistent sets are trivial (single 
valued) for regular connectivity. 
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Appendix B.3. Generation and stability of smaller set representations 

The problem with the preceeding method for generating a self-consistent basis for the 
1RSB method is twofold. Firstly we must begin by manipulating a potentially large set, 
the set Y bound , considering all possible states of C — 1 descendant fields. The method 
becomes unfeasible if C max is too large, greater than about five, since the number of 
field recombinations to be processed grows as a large power of C in (IB. 131) . For fields 
restricted to the two-state model description the scaling is more reasonable and creates 
no problem numerically up to C max ~ 10 in a brute force evaluation. However, it may 
be that the largest self-consistent set, which is the one we establish, may contain many 
fields irrelevant both to the cavity method solution and its local stability properties. 
There may exist some smaller yet sufficient sets that can be the basis of analysis. 

A sampling approach might be developed to address these problems. We identify a 
procedure that leads to a hierarchy of sets of increasing size, allowing greater control in 
the computational complexity. The method proposed relies on the recursive expansion 
of two sets. The first is T$ which is to be developed as the support for the solution, the 
second is a disjoint set of support perturbations Tp, which will be contracted to remove 
fields irrelevant to any local stability consideration. We define a mapping from these 
two sets to a third set T N as follows, 

r^v = T r (r 5 , r P ) = r 5 u {h E \h E = f({h E }) , h E _ x e r>, hf . . . h E _ 2 e r s } . (B.20) 

In the mapping only one of the descendants takes a field from the set of perturbations 
Tp. The set of fields in the perturbation set are considered bugs in an analogous manner 
to the type II stability analysis standard in 1RSB ^j, and do not interact. If under such 
a mapping T s = T N we say the set T s is linearly stable towards support bugs in the set 
T P . 

The procedure is therefore to select a simple set, using (IB. 141) to derive a self- 
consistent set Tg', then using (IB. 201) to discover the subset of perturbations that can 
survive the linearized recursion, Tp>. In principle we can now solve the 1RSB equations 
on the set Ts> and test the stability in the expanded space T N = T s > U Tpr. 

Assuming we find the solution is locally unstable, we can seek a larger support that 
might be sufficient to support the stable solution. Thus we can propose to repeat the 
process choosing Ts = T^. We summarize our method: 

1. Begin with an estimated set T$- 

2. Set r(* = °) = r^; derive a self consistent set recursively, applying ( 1B.14|) until 



convergence. Call the derived set T$'- If it does not converge repeat step 1 with a 
different set. 

3. Take a set T P of support bugs. 

4. Set r(*=°) = T P ; derive new sets recursively r(' +1 ) = T T , 2 (T S >,T^) \ T s > (lB~20|) . 
until convergence. Call the derived set Tpt. If it does not converge repeat step 3 
with a different set. 
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5. Return the pair {IV, T P /}, solve the energetic 1RSB method on T s > and test 1RSB 
instability across the support IV U T P i. 

6. If Tpi is the empty set, the support T s > is linearly stable against support bugs, 
terminate. Otherwise set Ts = IV U IV, and begin from step 1, to generate a new 
pair. 

Two types of stability analysis may be considered for the 1RSB solution, called type I 
and type II pi], and our method is able to identify the additional support necessary 
to rule out these linear instabilities that might go beyond the basis describing the 
solution. The art in the method comes from selecting sets that converge, avoiding 
cyclical behaviour that arises from symmetry breaking or some other effects. In fact, we 
find that the recursive generation procedure always produce converging sets provided 
one makes reasonable initial choices for the sets. 

The choice for IV is more involved, the simplest choice might be a single 
field meeting symmetry requirements and bounds. For ensembles with all nodes of 
connectivity at least |A| + 1 it is found that T s = {(0,0,0)} is self-consistent, and 
meets the requirements of stability against any perturbative set. This observation on 
the support alone implies that there must exist a locally stable trivial solution to the 
energetic RSB method in these ensembles. 

To generate non-trivial and asymmetric sets, those capable of describing a spin 
glass solution, we can consider the choice T basis. This can describe the 

symmetric solution, and with a reasonable choice of perturbations can be expanded to 
include the most relevant symmetry breaking fields. An alternative choice is Y^ rozen = 
{!// ,-J) .;■!>! .;)(// .;../} ,-!>j n) = (X h X 2 ,X 3 ),X i = ±X,0} with X -> oo, which 
is referred to as a frozen set since the effect of the infinite fields is to freeze the spin 
variables. This set initially violates the bounds, but converges upon iteration to a subset 
within the bounds. 

A choice for Tp, sufficient to consider any instability, is one that includes all possible 
self consistent fields, i.e. T 4 if it is known, or else Y bound , or all fields up to a symmetry 
constraint, i.e. T 2 or Y bound ' 2 . These are relatively large sets of perturbations, but since 
we need to consider only one of C — 1 descendants taking values from this large set, the 
combinatorial complexity is only linear in this number. 

Under the recursive set expansion procedure with either the frozen or asymmetric 
initial condition the sets converge at every stage, and following several iterations the 
process halts, in many cases with T s > = T 4 , or the trivial set IV = {(0,0,0)}. 
Interestingly from the frozen set, we find some other symmetry broken self-consistent 
sets, and also self-consistent sets smaller than T 4 for which the support is locally stable. 

Some results for the sizes of sets derived by our method are shown in the table, we 
denote as the self-consistent set derived from Y^ rozen using ( IB. 141) . and X + the set 
IV found by recursively applying our expansion method from an initial pair Ts = X 
and Tp = Y bound . After several iterations, with respect to the same perturbation set, 
a set of maximum size is found. A comparison of the support generated by various 
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Table Bl. The support for the energetic helds. 



methods is given in table IB1I Where the expanded sets are of equal size to T 4 they are 
the same set, where the set is of size one it is the trivial set {(0, 0, 0)}. We find that T 2+ 
always converges to one of these cases for A G [—5, 5] and C < 5, the trivial set only 
occurs for regular graphs in the special ratio (unlocked) regime A = C — 1. It is a trivial 
consequence, that the paramagnetic solution is locally stable in these cases at the level 
of energetic fields. However, since |T 4 | > 1 a more complicated solution may exist but 
cannot be found by a permutation procedure on the symmetric set of fields. 

Appendix B.3.1. Low energy simplifications of the field description In so far as the 
1RSB analysis is concerned, alternative simplifications can be proposed in order to 
allow manipulation of smaller sets. Firstly we can note that a very important special 
case of the 1RSB analysis involves the configuration parameter \i — > oo. The g-spin 
ancestor takes one of 4 states, and it is sufficient to consider only whether the state is 
excited (forbidden) or not excited, since any excited state contributes a vanishing weight 
to the pure-state distribution. The mappings of energetic fields (136|) can be modified 
and the notions of self-consistency and perturbative sets can be carried over to a set of 
size only 2 4 — 1. Since the size of this set does not scale with connectivity, or A, it is 
very transferrable between ensembles. 

More generally /i may be non-zero in the thermodynamically relevant solutions, 
and also we wish to describe metastable states. However, excited states are strongly 
penalized in general, and it can be expected that the phase space structure is mostly 
sensitive to only the local ground state structure. For this reason we can propose the 
auxiliary Hamiltonian (J6j). This removes all structures in the excitations, under the 
assumption that second level excitations and higher ones will have a weak effect on the 
energetic properties of the solution. 

For the auxiliary Hamiltonian we can carry out the same procedures as for the 
Quadratic Hamiltonian with significantly reduced set sizes. For this Hamiltonian 
a tighter set of bounds are immediately apparent, G {0, ±1}, and hence any 

symmetric solution has a support of three, and this can be combined with a modified 
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bound H b = C max + 1 (IB. 12[) which describes the symmetry broken solutions. A 
comparison of the set generation methods for the auxiliary Hamiltonian is shown in 
table IB~2l It can be seen that for all linear ensembles the largest self-consistent set (T 4 ), 
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Table B2. The support for the energetic helds, using the auxiliary Hamiltonian. 



used in the results section, is significantly smaller with the Auxiliary Hamiltonian. A 
comparison of results obtained based on the two Hamiltonians is made in section 17.41 

We considered also the method of this Appendix applied for ensembles with 
connectivity distribution supported on [2,C max ], for C max > 3, but not too large. 
This would describe a cut-Poisson ensemble with a truncated maximum connectivity, 
phenomena were found to be similar to the linear ensemble. 



Appendix C. Implementation of the Cavity Method 

Appendix C.l. Population dynamics initial condition for RS and 1RSB 

An important issue is the initial condition. For our finite temperature RS analysis we 
take fields to be 0(1) random numbers, with symmetry broken. 



The insight on the support gained by the methods of Appendix B can be applied 
equally to zero temperature RS and 1RSB field recursions. However, our zero- 
temperature RS results are developed with a population dynamics method, initializing 
every energetic field components independently from {0, ...,X}. The observables, 
stability and thresholds were not found to be particularly sensitive to the choice of 
X, provided a sufficient fraction of fields were distinct from (0,0,0), X is taken to be 
very large to obtain the results of this paper. A population of 10 4 fields is used in all 
RS results. 

In the 1RSB population dynamics, we consider a population of field distributions 
rather than a population of fields. In our experiments we chose 10 3 field distributions, the 



support of each distribution being the set T 4 identified in Appendix B Each member 
of our population is a vector of dimension |T 4 |, each component of the vector being 
the probability of the corresponding field hi. The components each distribution are 
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initialized by sampling uniformly from the interval [0,1], finally each distribution is 
normalized. 



Appendix C.2. Population dynamics convergence 

Beginning from the initial population we assume an ergodic regime will finally be reached 
in which the population moments are stable. For the replica symmetric solutions a 
minimum of 200 and maximum of 1200 iterations are allowed for convergence of the 
population. The spin glass order parameter qp in (T26]) and dimer magnetization run 
in ( 1271) . with sums replaced by population sampling, are expected to have a systematic 
trend during the transient phases of the population dynamics. This trend is towards 
in if the symmetric solution is the unique solution. We can consider the evolution of 
these observables over ~ 50 population updates, and allow additional iterations unless 
they have settled on a fixed mean up to statistical fluctuations. 

Similar criteria are applied to the 1RSB populations, the expressions of qp and mo 
are obtained by differentiating the 1RSB free energy, resulting in reweighted versions of 
( )26l) and ( )27l) under the hypothesis of many pure states. 



Appendix C. 3. The RS stability analysis for regular graphs 

Section 16.21 discusses a population dynamics approach to establishing solution stability, 
on a regular graph we have a local homogeneity that allows a simplification of the 
analysis. In this case the solution consistent with a single pure state, the RS 
or paramagnetic solution, is concentrated on h*. The joint distribution including 
perturbations can be described by 

P{h,5h) = 5(h-[h* + 8h])P{8h) (C.l) 

Suppose the distribution of perturbations is uniform over the population, that is, 
the shifts from (fjl^, fZ+i, ^j-n) — (0, 0, 0) are variations in the mean without any 
random component. We can then expand the mapping, a special case of (J34l) in which 
the perturbations on each kind of component are identical, and determine the new 
perturbation under iteration 

3 

5h{ +l = Y J M xy 5h t y (C.2) 

y=l 

where 

M xy = Jim dt{ f h}) (C.3) 

h-^h* dh y 

Under this mapping the perturbations may grow or decay exponentially, as determined 
by the largest real part of any eigenvalue associated with the 3 by 3 matrix M. 

It is sufficient to test only the stability of the mean and variance of the distribution 
to determine whether it is stable. Assuming that the mean of the distribution is stable 
we describe the remaining perturbations on ancestors and descendants by Gaussian 
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distributions parameterized by S. Each component T, yz describes the expectation value 
of the random quantity 5h y 5h z for the corresponding descendant or ancestor. Assuming 
the randomness to be of a homogeneous type on descendants, we can find the moments 
of the distribution f l4"5j) . and hence the mapping from descendant to ancestor covariance 

3 3 1 

= E E ^Ti^^SU (C4) 

J/=l 2 = 1 

As in the linear stability case, the principal eigenvalue (in the real part) of a 9 by 9 
matrix determines this non-linear stability, towards random zero mean perturbations. 



Appendix C.4- The 1RSB stability analyses overview 

The stability analysis is limited in our study to the ensembles C — 3, A = 1 with the 
quadratic and auxiliary Hamiltonians, and C = 4, A = 2 with the auxiliary Hamiltonian. 
These are locked problems of regular connectivity, and we consider the stability of the 



solution found with respect to the support T described in Appendix B Since for a 
regular graph ensemble there is no variation in the locally tree like neighborhoods, the 
1RSB solution can be described by a single distribution rather than a population of 
distributions. The standard type I and type II analyses are conducted [I]. We test the 
type I stability by adding perturbations to the distribution that evolve according to 
a linearized mapping. We monitor the first and second moments of the perturbations 
under iteration and find that the moments always decay under iteration, indicating 
stability. 

For type II instability, also called bug proliferation, we monitor how bugs in the 
distribution evolve on iteration. The bugs are described by infinitesimal probabilities 
P(h a ,h b ) that h a is "mistaken" for h^. Hence, on iteration the evolution of bugs is 
determined by a matrix |T 4 | 2 by |T 4 | 2 , relating the probabilities P(h a ,hb) before and 
after an iteration at a linearized level. The largest eigenvalue of this matrix determines 
the local stability. Rather than calculating this directly, we establish instability by 
randomly initializing a vector of the bug probabilities and measuring the sum of the 
variances of parturbations at successive iterations. Divergence of this sum implies at 
least one eigenvalue exceeds one, and the solution is unstable. We find that the sum 
diverges for every configurational parameter /i, and for both the quadratic and auxiliary 
Hamiltonian. Therefore the ground state solution, and metastable solutions, are locally 
unstable. Thus our solutions are type I stable, but type II unstable, for C = 3, A = 1 
and C = 4, A = 2. 



Appendix D. High Temperature Expansion in the A> 1 Limit 

This Appendix considers the case in which the next nearest neighbor interactions are 
much stronger than the nearest neighbor ones, thereby deriving the asymptote ( |65i) in 
the high temperature limit. Substituting the cavity probability distribution in (fl2"l) into 
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the recursion relation (J2J), we obtain the recursions identical to Eqns. ( 157)) to (159)) . with 

the energy D(Sj, Si, {h k ^.j}) substituted by the free energy F(Sj, Si, {h k ^j}) where 

2 



(d.i; 



F(Sj, S h {h k ^j}) = -ilogTr {5fc} exp \ (Si + Yfk=i S * 

+ E^i 1 ((A - J**) Sj - s k ] } . 

We note in passing that in the zero temperature limit, (ID.ip reduces to ( HU)) . Here we 
utilize Eq ( ID.lj) in the opposite limit of high temperature, in which the quadratic term 
in the Boltzmann factor is approximated as 



exp < — ft 



k=l 



c-i 



k=l 



This enables us to consider the effective Hamiltonian as one with purely nearest neighbor 
interactions, and the next nearest neighbors terms in (ID. 2}) as a thermodynamic average 
in the presence of such an effective Hamiltonian, that is, 

F(S 3 , Si, { W) « - \ logTr 5fc exp {-/? £ ((A - J*-i) Sj - /»*_►;) S k 



k=l 



(D.3) 



where the thermodynamic average in the last term is taken in the nearest neighbor 
model. Hence, apart from constant terms, we obtain 

1 

F{Sj, Si, {h k ^j}) « - - log cosh (/3 (A - J k ^) Sj - 



k=l 



£7-1 



5, ^ tanh (/3 (A - J fc ^) S 1 , - 
fc=i 

1 c_1 



(D.4) 

b 



fc<«=l 

x tanh (/3 (A- J^j) S s - Ptf^)] . 

To consider the stability of the paramagnetic phase, we consider cases that h^ i ,h h - 
are small. Following Eqs. ( 1571) to ( 1591) . we have 

£7-1 

-A + ^tanh(/3(A- J fc ^)) , (D.5) 



J 



fe=l 
£7-1 



-/3^/iLiSech 2 (/3(A- J*-,-)) 



(D.6) 



fc=i 



c-i 



£7-1 



*5-i = E h Li ~ E *U tanh(/3(A - J k ^)) 



k=l 



k=l 
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c-i 

+ /3^/iL,sech 2 (/3(A- J fc _>,))^tanh(/3(A- J k ^) . (D.7) 

k=l l^k 

Eq flD.5j) yields the high temperature paramagnetic state given by (Jj-^i, hj_^, hZ^) 
= (Jo, 0,0) where Jo = —A + (C — 1) tanh(2/3A). We are now ready to study the 
ferromagnetic and spin glass instabilities as described in Appendix C.3 The 3x3 
matrix ( 1C.3I) can be expanded to linear order in /3, to give the matrix M = M° + 0M 1 . 
Abbreviating t = tanh(2/3A) the two non-zero components of the matrix M° are 

M 3 ° 2 = (C-1) and M$ 3 = -(C-l)t, (D.8) 

implying that the non-zero eigenvalue at the leading order is eo = — (C — l)t. This 
eigenvalue is perturbed according to the elements of M 1 , whose non-zero elements are 

M\ x = M\ 3 = -{G - 1)(1 - t 2 ) and = (C - 1)(2C - 3)t(l - t 2 ) . (D.9) 

To establish the perturbation in the eigenvalue, we consider the matrix M' = M — e^I, 
I being the identity matrix. The relevant eigenvalue of M' is of the order /3. Hence it 
suffices to consider the characteristic polynomial of M' to the first order in /3, given by 

C(x) = (-e ) 2 -x)- (-eo)/3M 2 1 3 M 3 2 . (D.10) 

The eigenvalue is 

/9ei = -(C - 1)(1 - t 2 )[l + (2C - 3)t 2 } . (D.ll) 

The ferromagnetic instability is determined by equating the eigenvalue eo + 0e\ to 1. 
Expanding A about its zeroth order value, and inverting its dependence on temperature, 
we find 

2 / C 2 — 2 \ 

Tc = "atanh((C-l)-i) ( A + 2(C^T) ) ' (D ' 12) 
We find that the Curie temperature is always reduced by the presence of next nearest 
neighbor antiferromagnetic interactions, and by 0(C 2 ). 

To calculate the spin glass instability, we determine the largest eigenvalue of the 
9x9 matrix M^^i) = MikMjt/(C — 1), as was outlined in Appendix C.3 The leading 
order of the eigenvalue is determined to be eo = M® 33 ^ , 33 s = (C — l)t 2 . To 0(/3), the 
eigenvalue e + /3ei is determined by the equation 

(-e ) 8 (C - 1)- 1 (2/3M 3 3M 3 1 3 - e x ) - (-e ) 7 (C - 1 ) _ 1 2 /3 M 32 ( M 3 3 ) 2 M\ 3 = , (D.13) 
whose solution is 

e x = 2(C - 1)(1 - t 2 )[l + (2C - 3)t 2 ] . (D.14) 

By equation the eigenvalue to 1, the critical curve ( 1651) is thereby obtained. 
Acknowledgements 



NNNI models on random graphs 



54 



Acknowledgments 

We thank Bill Yeung and David Saad for fruitful discussions, and Stefan Boettcher for 
providing source code that was adapted for this investigation. The work is supported 
by Research Grants Council of Hong Kong (grant numbers 604008 and 605010). 



References 



M. Mezard, G. Parisi, and M.A Virasoro. Spin Glass Theory and Beyond. World Scientific, 
Singapore, 1987. 

H. Nishimori. Statistical Physics of Spin Glasses and Information Processing. Oxford Science 

Publications, Oxford, UK, 2001. 
F. Krzakala, A. Montanari, F. Ricci-Tersenghi, G. Semerjian, and L. Zdeborova. Gibbs states and 

the set of solutions of random constraint satisfaction problems. PNAS, 104(25):10318-10323, 

2007. 

A. Montanari and F. Ricci-Tersenghi. On the nature of the low-temperature phase in discontinuous 

mean-field spin glasses. Eur. Phys. J. B, 33:339, 2003. 
W. Seiko. The ANNNI model- theoretical analysis and experimental application. Phys. Rep., 

170(4):213-264, 1988. 

C. Domb and R.B. Potts. A two-dimensional model with first and second interactions. Proc. Roy. 

Soc. Lond. A., 210:125 - 141, 1951. 
J. Stephenson and D. D. Betts. Ising model with antiferromagnetic next-nearest-ncighbor coupling. 

ii. ground states and phase diagrams. Phys. Rev. B, 2(7):2702-2706, 1970. 
J.N. Reimers and A.J. Bcrlinsky. Order by disorder in the classical Heisenberg kagome 

antifcrromagnct. Phys. Rev. B, 48(13) :9539-9554, Oct 1993. 
M.H. Jensen and P. Bak. Mean-field theory of the three-dimensional anisotropic Ising model as a 

four-dimensional mapping. Phys. Rev. B, 27(ll):6853-6868, Jun 1983. 
J. Vannimcnus. Phase diagram of an Ising model with competitive interactions on a Husimi tree 

and its disordered counterpart. Z. Phys. B, 43:141 - 148, 1981. 
J. G. Moreira and S. R. Salinas. Modulated structures in the Ising model with competing 

interactions on the Cayley tree. Phys. Rev. B, 47:778-786, January 1993. 
C. S. O. Yokoi, M. J. de Oliveira, and S. R. Salinas. Strange attractor in the Ising model with 

competing interactions on the Cayley tree. Phys. Rev. Lett., 54(3):163-166, Jan 1985. 
S. Inawashiro, C.J. Thompson, and G. Honda. Ising model with competing interactions on a 

Cayley tree. J. Stat. Phys., 33(2):419 - 436, 1983. 
S. Katsura and M. Takizawa. Bethe lattice and the Bethe approximation. Prog. Theor. Phys., 

51(l):82-98, 1974. 

N.N. Ganikhodjaev, C.H. Pah, and M.R.B. Wahiddin. Exact solution of an Ising model with 

competing interactions on a Cayley tree. J. Phys. A, 36(15):4283-4289, 2003. 
C.R. da Silva and S. Coutinho. Ising model on the Bethe lattice with competing interactions up 

to the third-nearest-neighbor generation. Phys. Rev. B, 34(ll):7975-7985, Dec 1986. 
N. Ganikhodjaev, F. Mukhamedov, and C.H Pah. Phase diagram of the three states potts model 

with next nearest neighbour interactions on the Bethe lattice. Phys. Lett. A, 373(1):33 - 38, 

2008. 

R. Mclin and S. Pcysson. Spin glass behavior upon diluting frustrated magnets and spin liquids: 

a Bethe-Peierls treatment. Eur. Phys. J. B, 14(1):169-176, 2000. 
J.L. Monroe. Phase diagrams of Ising models on Husimi trees ii. pair wand multisite interaction 

systems. J. Stat. Phys., 67(5):1185-1200, 1992. 
P Chandra and B Doucot. Spin liquids on the Husimi cactus. Journal of Physics A: Mathematical 

and General, 27(5): 1541, 1994. 



NNNI models on random graphs 



55 



[21] M. Ostilli, F. Mukhamedov, and J.F.F. Mendes. Phase diagram of an Ising model with competitive 

interactions on a Husimi tree and its disordered counterpart. Physica A: Statistical Mechanics 

and its Applications, 387(12):2777 - 2792, 2008. 
[22] J.G. Moreira and S.R. Salinas. Ising model with third-neighbour interactions on the Cayley tree. 

J. Phys. A, 20(6):1621, 1987. 
[23] K.Y.M. Wong and D. Saad. Minimizing unsatisfaction in colourful neighbourhoods. J. Phys. A, 

41(32):324023 (25pp), 2008. 
[24] S. Bounkong, J. van Mourik, and D. Saad. Coloring random graphs and maximizing local diversity. 

Phys. Rev. E, 74(5):057101, 2006. 
[25] A. Pelizzola, M. Pretti, and J. van Mourik. Palette-colouring: a belief-propagation approach. 

arXiv:1104.4024, 2011. 

[26] M. Kearns, M. Littman, and S. Singh. Graphical models for game theory. In Proceedings of the 

Seventeenth Conference Annual Conference on Uncertainty in Artificial Intelligence (UAI-01), 

pages 253-260, San Francisco, CA, 2001. Morgan Kaufmann. 
[27] Ramezanpour, A., Realpe-Gomez, J., and Zecchina, R. Statistical physics approach to graphical 

games: local and global interactions. Eur. Phys. J. B, 81(3):327-339, 2011. 
[28] M. Mezard and G. Parish The cavity method at zero temperature. J. Stat. Phys., lll(l-2):l-34, 

2003. 

[29] J. S. Yedidia, W. T. Freeman, and Y. Weiss. Constructing free energy approximations and 

generalised belief propagation algorithms. Technical Report TR2002-35, Mitsubishi Electric 

Research Laboratories, 2002. 
[30] F.R. Kschischang, B.J. Frey, and Hans-Andrea Loeliger. Factor graphs and the sum-product 

algorithm. IEEE Trans, on Info. Theory, 47(2):498-518, 2001. 
[31] L. Zdcborova and M. Mezard. Constraint satisfaction problems with isolated solutions are hard. 

J. Stat. Mech., 2008(12):P12004, 2008. 
[32] A. Janson, T. Luczak, and A. Rucinski. Random Graphs. John Wiley & sons, New York, NY, 

USA, 2000. 

[33] S. Boettcher and A. Pcrcus. Nature's way of optimizing. Artif. Intell, 119(l-2):275-286, 2000. 
[34] S. Boettcher. Numerical results for ground states of mean-field spin glasses at low connectivities. 

Phys. Rev. B, 67:060403, 2003. 
[35] R.J. Baxter. Exactly Solved Models in Statistical Mechanics. Academic Press, New York, NY, 

USA, 1982. 

[36] M. Mezard and G. Parisi. The Bethe lattice spin glass revisited. Eur. Phys. J. B, 20(2):217-233, 
2001. 

[37] A. Montanari, G. Parisi, and F. Ricci-Tersenghi. Instability of one-step replica-symmetry-brokcn 
phase in satisfiability problems. Journal of Physics A: Mathematical and General, 37(6):2073, 
2004. 

[38] M. Bayati, C. Borgs, J. Chayes, and R. Zecchina. On the exactness of the cavity method for 
weighted b-matchings on arbitrary graphs and its relation to linear programs. Journal of 
Statistical Mechanics: Theory and Experiment, 2008(06):L06001, 2008. 

[39] L. Zdeborova and M. Mezard. The number of matchings in random graphs. Journal of Statistical 
Mechanics: Theory and Experiment, 2006(05):P05003, 2006. 

[40] I. Kanter and H. Sompolinsky. Mean-field theory of spin-glasses with finite coordination-number. 
Phys. Rev. Lett, 58(2):164-167, 1987. 



